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ABSTRACT 

We derive a new equation of state (EoS) for neutron stars (NS) from the outer crust to the core based on modern microscopic 
calculations using the Argonne Uig potential plus three-body forces computed with the Urbana model. To deal with the inhomogeneous 
structures of matter in the NS crust, we use a recent nuclear energy density functional that is directly based on the same microscopic 
calculations, and which is able to reproduce the ground-state properties of nuclei along the periodic table. The EoS of the outer crust 
requires the masses of neutron-rich nuclei, which are obtained through Hartree-Fock-Bogoliubov calculations with the new functional 
when they are unknown experimentally. To compute the inner crust, Thomas-Fermi calculations in Wigner-Seitz cells are performed 
with the same functional. Existence of nuclear pasta is predicted in a range of average baryon densities between =^0.067 fm“^ and 
-0.0825 fm“^, where the transition to the core takes place. The NS core is computed from the new nuclear EoS assuming non-exotic 
constituents (core of npep matter). In each region of the star, we discuss the comparison of the new EoS with previous EoSs for the 
complete NS structure, widely used in astrophysical calculations. The new microscopically derived EoS fulfills at the same time a NS 
maximum mass of 2 with a radius of 10 km, and a 1.5 Mq NS with a radius of 11.6 km. 
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1. Introduction 

Neutron stars (NS) harbor unique conditions and phenomena 
that challenge the physical theories of matter. Beneath a thin 
stellar atmosphere, a NS interior consists of three main re¬ 
gions, namely, an outer crust, an i nner crust, and a core, each 
one featuring a different physics ( Shapiro & Teukolskvl 1 19831: 
iHaensel et all 1200^ IChaniel & Haensell 20081) . The core is the 
internal region at densities larger than 1.5x10'"^ g/cm^, where 
matter forms a homogeneous liquid composed of neutrons plus 
a certain fraction of protons, electrons, and muons that maintain 
the system in p equilibrium. Deep in the core, at still higher den- 
sities, strange baryons and even deconfined quar ks may appear 
(IShaniro & Teukolskvil 19831 IHaensel et alJl2()()^ . Moving from 
the core to the exterior, density and pressure decrease. When the 
density becomes lower than approximately 1.5x10^^^ g/cm^, mat¬ 
ter inhomogeneities set in. The positive charges concentrate in 
individual clusters of charge Z and form a solid lattice to min¬ 
imize the Coulomb repulsion among them. The lattice is em¬ 
bedded in a gas of free neutrons and a background of electrons 
such that the whole system is charge neutral. This region of the 
star is called the inner crust, where the nuclear structures may 
adopt non-spherical shapes (generically re ferred to as nuclear 
pasta) in order to minimize their energy ( Bavm et al. jjTlal : 
iRavenhall et alJfT983l iLorenz et al.lll99!M lOvamatsul 1993|) . At 
lower densities, neutrons are finally confined within the nuclear 
clusters and matter is made of a lattice of neutron-rich nuclei 
permeated by a d egenerate electron g as. This region is known 
as the outer crust (iBavm et alJll971a) and extends, from inside 
to outside, from a neutron drip density of about 4xl0" g/cm^ 
to a density of about 10"^ g/cm^. Most of the mass and size of 


a NS are accounted for by its core. Although it is only a small 
fraction of the star mass and radius, the crust plays an impor¬ 
tant role in various observed astrophysical phenomena such as 
pulsar glitches, quasiperiodic oscillations in soft gamma-ray re¬ 
peaters (SGR), and thermal relaxation in soft X-ray transients 


(SXT) (IHaensel et al.l2006l:IChamel & HaenseJ200alPirol2005 


Sttohma^er_&^^tts|j200d Steiner&^^^^itts[ 2009t ISotani et al.l 


20121 : iNewton et alJ l2013t iPiekarewicz et al. l201'4 . which de¬ 
pend on the departure of the star from the picture of a homoge¬ 
neous fluid. Recent studies suggest that the existence of nuclear 
pasta layers in the NS crust may limit the rotational speed of 
pulsars and may be a p ossible origin of the lack of X-ray pulsa rs 
with long spin periods (iPons et al]l2013l [Horowitz et alJl2015l) . 

The equation of state (EoS) of neutron-rich matter is a basic 
input needed to compute most properties of NSs. A large body of 
experimental data on nuclei, heavy ion collisions, and astrophys¬ 
ical observations has been gathered and used along the years to 
constrain the nuclear EoS and to understand the structure and 
properties of NS. Unfortunately, a direct link of measurements 
and observations with the underlying EoS is very difficult and a 
proper interpretation of the data necessarily needs some theoret¬ 
ical inputs. To reduce the uncertainty on these types of analyses, 
it is helpful to develop a microscopic theory of nuclear matter 
based on a sound many-body scheme and well-controlled basic 
interactions among nucleons. To this end, it is of particular rele¬ 
vance to have a unified theory able to describe on a microscopic 
level the complete structure of NS from the outer crust to the 
core. 

There are just a few EoSs devised and used to describe 
the whole NS within a unified framework. It is usually as¬ 
sumed that the NS crust has the structure of a regular lattice 
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that is treated in the Wigner-Seitz (WS) approximation. A par¬ 
tially phenomen ological approach was de veloped by Lattimer 
and Swesty (LS) (iLattimer & Swest^ll991h . The inner crust was 
computed using the compressible liquid drop model (CLD M) 
introduced by Baym, Bethe, and Pethick (iBavm et al.lll97l3) to 
take into account the effect of the dripped neutrons. In the LS 
model the EoS is derived from a Skyrme nu clear effective force. 
There are different v ersions of this EOS (ILattimer & Swest^ 
1199 It lLattimerll2015h . each one having a different i ncompress¬ 
ibility. A nother EOS was developed by Shen et al. (IShen et al.l 
ll998bllaL ISumivoshill20 1 5l) based on a nuclear relativistic mean 
field (RME) model. The crust was described in the Thomas- 
Eermi (TE) scheme using the variational method with trial pro¬ 
files for the nucleon densities. The LS and Shen EoSs are widely 
used in astrophysical calculations for both neutron stars and su¬ 
pernova simulations due to their numerical simplicity and the 
large range of tabulated densities and temperatures. 

Douchin and Haensel (DH) (iDouchin & Haens'dll2001l) for¬ 
mulated a unified EoS for NS on the basis of the SLy4 Skyrme 
nuclear effective force (IChabanat et al.lfT998h . where some pa¬ 
rameters of the Skyrme interaction were adjuste d to reproduce 
the W iringa et al. calculation of neutron matter (IWiringa et al.l 
Il988h above saturation density. Hence, the DH EoS contains 
certain microscopic input. In the DH model the inner crust 
was treated in the CLDM approach. More recently, unified 
EoSs for NS have bee n derived by the Bru s sels-Montreal group 
iPearson et al.l 1201 2t lEantina et al.l 120131: 
They ar e based on the BS k family of 


(Charnel et al.l 

201 1[ 

IPotekhin et al.l 

2013b 


Skyrme nuclear effectiye forces (iGoriely et al.l2010l) . Each force 
is fitted to the known masses of nuclei and adjusted among 
other constraints to reproduce a different microscopic EoS of 
neutron matter with different stiffness at high density. The in¬ 
ner crust is treated in the extended Thomas-Eermi approach 
with trial nucleon density profiles including perturbatiyely shell 
corrections for protons yia the Strutinsky integral method. An¬ 
alytical fits of these neutron-star EoSs haye been constructed 
in or der to facilitate thei r inclusion in astrophysical simula¬ 
tions (iPotekhin et al.ll20T^ . Quantal Hartree calculat ions for the 
NS cru st haye been systematically performed by (IShen et al.l 
1201 IbO . This approach uses a yirial expansion at low den¬ 
sity and a RME effectiye interaction at intermediate and high 
densities, and the EoS of the whole NS has been tabulated for 
different RME parameter sets. Also recently, a complete EoS 
for sup ernpya matter has been deyeloped within the statistical 
model dHempel & Schaffner-Bielic h||2010|). We shall adopt here 


the E o S of the BSk21 mod el ([Charnel eta 


2012; lEantina et al.l 120131: IPotekhin et al 


201 It IPearson et al. 


m 

II 


2013t IGoriely et al. 


20101) as a representatiye example of contemporary EoS for 
the complete NS struc ture, and a comparison with the other 
EoSs of the BSk far nily ( Charnel et al.|201 U lPearson et al.l20i:^ 
lEantina et al.|| 201^ IPotekhin et al. 20131) and the RME family 
( Shen et al.ll201 Ibll^ shall be left for future study. 

Our aim in this paper is to obtain a unified EoS for neu¬ 
tron stars based on a microscopic many-body theory. The nu¬ 
clear EoS of the model is deriyed in the Brueckner-Hartree- 
Eock (BHE) approach from the bare nucleon-nuc leon (NN) in¬ 
teract ion in free space (Argonne wig potential (IWiringa et al.l 
1 19951) ) with inclusion of three-body forces (TBE) among nucle¬ 
ons. In the current state of the art of the Brueckner approach, 
the TBEs are reduced to a density-dependent two -body force 
by ay eraging oyer the third nucleon in the medium dBaldo et al.l 
Il997h . We employ TBEs based on the Urbana model, consist¬ 
ing of an attractiye term from two-pion exchange and a repul- 
siye phenomenological central term, to reproduce the satura¬ 


tion point (Schiayilla et al.lfT9^ iBaldo et ^Il997t iBaldo et ^ 
1201 21: iTaranto et al.ll201^ The corresponding nuclear EoS for 
symmetric and asymmetric nuclear matter fulfills seyeral re¬ 
quirements impos ed both by heayy io n collisions and astrophys¬ 
ical obseryations dTaranto et al.ll201^ . 


Recently the connection between two-body and three- 
body forces within the meson-nucleon theory of the nu¬ 
clear in teraction has been extensiyely discussed and deyel- 


oped in( 

Grange et al.ll989HZuo et al.l2002HLi & Schulzel2008l 

iLi et al.l 

20081). At present the theoretical status of microscop- 


ically deriyed TBEs is still incipient, howeyer a tentatiye ap¬ 
proach has been proposed using the same meson-exchange pa¬ 
rameters as the underlying NN potential. Results haye been 
obtained with the Argonne fig, the Bonn B, and t he Ni¬ 
jmegen 93 potentials dLi & Schulzell20d^ iLi et al.ll2008l) . More 


recently the chiral expansion theory t o the nucleon interac- 


tion has been extensiyely deyeloped (IWeinbergI 19681 19901 

1991 

11992; Entem & Machleidtl 120031 Valderrama & Phillios 

2015 

; Leutwylei 19941; Epelbaum et al.l 20091; Otsukaetal. 

2010 

; Holt et al.l2013t Hebeler et al.l201 It 

Drischler et al.l2014l; 

Hebeler & Schwenk 

20101 Carbone et al.l 

2013; Ekstrom et al. 

2013; Coraggio et al. 

2014h. This approach is based on a deeper 


leyel of the strong interaction theory, where the QCD chi¬ 
ral symmetry is explicitly exploited in a low-momentum ex¬ 
pansion of the multi-nucleon interaction processes. In this ap¬ 
proach multi-nucleon interactions arise naturally and a hierar¬ 
chy of the different orders can be established. Despite some 
ambig uity in the parametrization of the force dEkstrom et al.l 
I 2 OI 3 I) and some diffi culty in the treatment of many-body systems 
( Laehde et al.ll2013l) . the method has marked a great progress in 
the microscopic theory of nuclear system s. Indeed it turns out 
(ICoraggio et al.ll20T4tlEkstr6m et al.l(2015l) that within this class 
of interactions a compatible treatment of few-nucleon systems 
and nuclear matter is possible. Howeyer, this class of NN and 
TBE (or multi-body) interactions is deyised on the basis of a 
low-momentum expansion, where the momentum cut-off is fixed 
essentially by the mass of the p-meson. As such they cannot 
be used at density well aboye saturation, where we are also go¬ 
ing to test the proposed EoS. In any case, the strength of TBEs 
is model dependent. In particular, the role of TBEs appears to 
be marginal in the quark- meson model of the NN interaction 
(IB aldo & Eukuka\^l20 1 4l) . 


A many-body calculation of the inhomogeneous structures 
of a NS crust is currently out of reach at the leyel of the BHE 
approach that we can apply for the homogeneous matter of the 
core. In an attempt to maximize the use of the same microscopic 
theory for the description of the complete stellar structure, we 
employ in the crust calculations the recently deyeloped BCPM 
(Barce lona-Catania-Paris-Madrid) nucl ear energy density func¬ 
tional (IBaldo etal.ll200^ I 20 M BOTl . The BCPM functional 


has been obtained from the ab initio BHE calculations in nuclear 
matter within an approximation inspir ed by the Kohn-Sham for¬ 
mulation of density functional theory (iKohn & Shaml[l96^ . In¬ 
stead of starting from a certain effectiye interaction, the BCPM 
functional is built up with a bulk part obtained directly from the 
BHE results in symmetric and neutron matter yia the local den¬ 
sity approximation. It is supplemented by a phenomenological 
surface part, which is absent in nuclear matter, together with the 
Coulomb, spin-orbit, and pairing contributions. This energy den¬ 
sity functional constructed upon the BHE calculations has a re¬ 
duced set of four adjustable parameters in total and describes the 
ground-state properties of finit e nuclei similarly successfully a s 
the Skyrme and Gogny forces (IBaldo et al.ll2008dl2010il2013l) . 
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We model the NS crust in the WS approximation. To com¬ 
pute the outer crust, we take the masses of neutron-rich nuclei 
from experiment if they are measured, and perform deformed 
Hartree-Fock-Bogoliubov (HFB) calculations with the BCPM 
energy density functional when the masses are unknown. To 
describe the inner crust, we perform self-consistent Thomas- 
Fermi calculations with the BCPM functional in different peri¬ 
odic configurations (spheres, cylinders, slabs, cylindrical holes, 
and spherical bubbles). In this calculation of the inner crust, by 
construction of the BCPM functional, the low-density neutron 
gas and the bulk matter of the high-density nuclear structures are 
not only consistent with each other, but they both are given by the 
microscopic BHF calculation. We find that nuclear pasta shapes 
are the energetically most favourable configurations between a 
density of ^ 0.067 fm“^ (or 1.13x10'^ g/cm^) and the transi¬ 
tion to the core, though the energy differences with the spherical 
solution are small. The NS core is assumed to be composed of 
npejj. matter. We compute the EoS of the core using the nuclear 
EoS from the BHE calculations. In the past years, the BHE ap¬ 
proach has b een extended in order to in clude the hyperon degrees 
of freedom (iBaldo et al.lll998ll2000ah . which play an important 
role in the study of neutron-star matter. However, in this work we 
are mainly interested on the properties of the nucleonic EoS, and 
therefore we do not consider cores with hyperons or other exotic 
components. In each region of the star we critically compare the 
new NS EoS with the results from various semi-microscopic ap¬ 
proaches mentioned above, where the crust and the core were 
calculated within the same theoretical scheme. Lastly, we com¬ 
pute the mass-radius relation of neutron stars using the unified 
EoS for the crust and core that has been derived from the micro¬ 
scopic BHE calculations in nuclear matter. The predicted maxi¬ 
mum mass and radii are compatible with the recent astrophysical 
observations and analyses. We reported some prel iminary results 
about the new EoS recently in (iBaldo et al.ll2014t) . 

In Sec. |2] we summarize the microscopic nuclear input to 
our calculations and the derivation of the BCPM energy density 
functional for nuclei. We devote Sec. [3]to the calculations of the 
outer crust. In Sec.|4]we introduce the Thomas-Eermi formalism 
for the inner crust with the BCPM functional, and in Sec. |5]we 
discuss the results in the inner crust, including the pasta shapes. 
In Sec. 0 we describe the calculation of the EoS in the core and 
obtain the mass-radius relation of neutron stars. We conclude 
with a summary and outlook in See.]?] 


2. Microscopic input and energy density functionai 
for nuciei 

As the microscopic BHE EoS of nuclear matter underlies our 
formulation, we start by summarizing how the EoS of nuclear 
matter is obtained in the BHE method and how it has been used 
to construct a suitable energy density functional for the descrip¬ 
tion of finite nuclei. 

The nuclear EoS of the model is derived in the framework 
of the Brueckner-Bethe-Goldstone theory, which is based on a 
linked clust er expansion of the energy per nucleon of nuclear 
matter (see (lBaldolll999t) . chapter 1 and references therein). The 
basic ingredient in this many-body approach is the Brueckner 
reaction matrix G, which is the solution of the Bethe-Goldstone 
equation 


G[n; to] - V + 


Z 

ka,h 


V 


\^a^h} Q iha^b\ 

CO - eika) - e(kb) 


G[n; to ], 


( 1 ) 


where v is the bare NN interaction, n is the nucleon number den¬ 
sity, and to is the starting energy. The propagation of intermedi¬ 
ate baryon pairs is determined by the Pauli operator Q and the 
single-particle energy e{k), given by 


e{k) = e{k-,n) = — + U{k-,n). 
2m 


( 2 ) 


We note that we assume natural units ti - c - I throughout the 
paper. The BHE approximation for the single-particle potential 
t/(k; n) using the continuous choice is 

U(k; n)^Yj {kk'\G[n-, e{k) + e{k')]\kk')^ , (3) 

k'<kf 


where the matrix element is antisymmetrized, as indicated by 
the “a” subscript. Due to the occurrence of U{k) in Eq. (I2]i, the 
coupled system of equations O to Q must be solved in a self- 
consistent manner for several Eermi momenta of the particles 
involved. The corresponding BHE energy per nucleon is 

X {kAGMk) + e{k')]\kk\. (4) 

k,k'<kF 

In this scheme, the only input quantity we need is the bare NN in¬ 
teraction V in the Bethe-Goldstone equation ([1]). The nuclear EoS 
can be calculated with good accuracy in the Brueckner two-hole- 
line approximation with the continuous choice for the single¬ 
particle potential, since the results in this scheme are quite close 
to the cal culations which i nclude also the three-hole-line con¬ 
tribut ion dSong et al.lll998t IBaldo et ^l2000bt IBaldo & Burgiol 
l2001h. In the presen t work, we use the Argonne potential 
( Wiringa et al.l 19951) as the two-nucleon interaction. The depen¬ 
dence on the NN interaction, as well as a comparison with other 
many-body approa ches, has been systematically investigated in 
(IBaldo et al.ll2012h . 

One of the well-known results of several studies, which 
lasted for about half a century, is that non-relativistic calcula¬ 
tions, based on purely two-body interactions, fail to reproduce 
the correct saturation point of symmetric nuclear matter, and 
three-body forces among nucleons are needed to correct this 
deficiency. In our approach the TBE is reduced to a density- 
dependent two-body force by averaging over the position of 
the third particle, assuming that the probability of having two 
particles at a given distan ce is reduced acco rding to the two- 
body correlation function (IBaldo et al.lll997l) . In this work we 
use a phenomenological approach based on the so-called Ur- 
bana model, which consists of an attractive term due to two- 
pion exchange with excitation of an intermedi ate A resonance, 
and a repulsive phenomenological central term (ISchiavilla et al.l 
Il986h . Those TBEs produce a shift of the saturation point (the 
minimum) of about -i-l MeV in energy. This adjustment was 
obtained by tuning the two parameters contained in the TBEs, 
and was performed to g et an optimal saturation point (for details 
see dBaldo et al.l 1 1997^ 1. The calculated nuclear EoS conforms 
to several constraints required by the phenomeno logy of heavy 
ion co llisions and astrophysical observational data dTaranto et al.l 
1201 3h . 

Eor computational purposes, an educated polynomial fit is 
performed on top of the microscopic calculation of the nuclear 
EoS including a fine tuning of the two parameters of the three- 
body forces such that the saturation point be E/A = -16 MeV 
at a density no = 0.16 fm“^. The interpolating polynomials for 
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Table 1. Coefficients of the polynomial fits E/A of the EoS of symmet¬ 
ric matter and neutron matter, see Eq. (0. 


k 

au (MeV) 

bk (MeV) 

1 

-73.292028 

-34.972615 

2 

49.964912 

22.182307 

3 

-18.037608 

-7.151756 

4 

3.486179 

1.790874 

5 

-0.243552 

-0.169591 


symmetric nuclear matter and pure neutron matter are written as 


=(f 
=(f L, 



(5) 


where no„ = 0.155 fm“^. The values of the coefficients of the 
interpolating polynomials (BCP09 version (iBaldo et al.ll2010l) 
of the parameters) are given in Table [1] This fit is valid up to 
density around n = 0.4 fm“^, and it is used only for conve¬ 
nience. For higher density, as the ones occurring in NS cores, 
we use a direct numerical interpolation of the computed EoS, 
using the same procedur e and functional form as in reference 
dBurgio & Schulzjl2010h for a set of accurate BHF calculations. 
We found that the polynomial form of the EOS and the high 
density fit join smoothly in the interval 0.3 - 0.4 fm“^, for both 
the energy density and the pressure of the fi stable matter. The 
overall EOS so obtained will be reported in the corresponding 
Tables below and used for NS calculations. The properties of in¬ 
finite nuclear matter at saturation are collected in Table |2l and 
we see that their values agree very well with the known empir¬ 
ical values. The symmetry energy and the corresponding slope 
parameter L are two important quantities closely related to var¬ 
ious properties of neutron stars and to the thickne ss of the neu¬ 
tron skin of nuclei dHorowitz & Piekarewiczl2001h . It can be no¬ 
ticed that the predicted values for EsymCHo) and L lie within the 
recent constraints deri ved from the analys i s of different astro- 
physical observations dNewton_et^ 200 ISteiner et alJ 120101 : 
iLattimer & Lim 2013: Hebeler et al. 1201 3l) and nuclear experi- 


ments dChen et al.ll2010l : lTsang et alJl2012tlVinas et'ani201'4 . 


The BHF result can be directly employed for the calculations 
of the NS liquid core, where the nuclei have dissolved into their 
constituents, protons and neutrons. However, the description of 
finite nuclei and nuclear structures of a NS crust is not man¬ 
ageable on a fully microscopic level. Indeed, the only known 
tractable framework to solve the nuclear many-body problem 
in finite nuclei across the nuclear chart is provided by density 
functional theory. In order to describ e finite nuclei, the BCPM 
energ y density functional was built dBaldo et al.ll2008R 120101. 
12013h based on the same microscopic BHF calculations pre¬ 
sented before. The BCPM functional is obtaine d within an ap¬ 
proxi mation inspired by the Kohn-Sham method dKohn & ShamI 
Il965h . In this approach the energy is split into two parts, the 
first one being the uncorrelated kinetic energy, while the second 
one contains both the potential energy and the correlated part 
of the kinetic energy. An auxiliary set of A orthonormal orbitals 
1 , 0 ,(r, cr, q) is introduced (where A is the nucleon number and cr 
and q are spin and isospin indices), allowing one to write the 
one-body density as if it were obtained from a Slater determi¬ 
nant as n{r) - Yii,a-,q ?)l^, and the uncorrelated kinetic 
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energy as 

^ 0 =^^ ^\Vipi{r,a-,q)\^dr. (6) 


To deal with the unknown form of the correlated part of the 
energy density functional, a strategy often followed in atomic 
and molecular physics is to use accurate theoretical calculations 
performed in a uniform system which finally are parametrized 
in terms of the one-body density. We use a similar approach, 
and we apply it to the nuclear many-body problem to obtain the 
BCPM functional. For this purpose we split the interacting nu¬ 
clear part of the energy {E^) into bulk and surface contributions 
(i.e. En - E^^^ + F™'^). We obtain the bulk contribution 
directly from the ab initio BHF calculations of the uniform nu¬ 
clear matter system by a local density approximation. Namely, 
in our approach depends locally on the nucleon density 
11 - n„ + lip and the asymmetry parameter/? = {n„-np)/(n„ + Up) 
and reads as 

F“'‘K,n„] = J [P.(n)(l-/?2) + P„(n)yS2]nr/r, (7) 


where the polynomials Pj(n) and P„(n) in powers of the one- 
body density have been introduced in Eq. (I5]l. 

In addition to the bulk part, also surface. Coulomb, and spin- 
orbit contributions which are absent in nuclear matter are neces¬ 
sary in the interacting functional to properly describe finite nu¬ 
clei. We make the simplest possible choice for the surface part 
of the functional by adopting the form 

Fr^[np,«n] = «9(r)%?'(r - r')?v(r')£/rc/r' 

- J'/ig(r)ng^(r)dr Jug^g^ir'/dr', ( 8 ) 

where q = n,p for neutrons and protons. The second term in ® 
is subtracted to not contaminate the bulk part determined from 
the microscopic nuclear matter calculation. For the finite range 
form factors we take a Gaussian shape Vg^g'(r) - Vg^g>e~'^^°‘\ 
with three adjustable parameters: the range a, the strength Vp^p = 
y„ „ = Vl for like nucleons, and the strength V„g, - Vp_„ = Vu 
for unlike nucleons. The Coulomb contribution to the functional 
is the sum of a direct term and a Slater exchange term computed 
from the proton density: 



np(r)np(r') 

|r-r'| 


3e^ 

dr dr' -— 

4 



(9) 


As in Skyrme and Gogny forces, the spin-orbit term is a zero- 
range interaction o = iWo(<T, -H CTj) x [k' x b(r, - r^jk], whose 
contribution to the energy reads 




[«(r)VJ(r) + «„(r)VJ„(r) + n„(r)VJ„(r)]ufr, (10) 


where the uncorrelated spin density 3g for each kind of nucleon 
is obtained using the auxiliary set of orbitals as 


J?(r) = -i ^ ^f)*(r,cr,^)[V(,c,(r,cr',^)x<cr|cr|cr')]. (11) 

iy(T,cr' 


Thus, the total energy of a finite nucleus is expressed as 
£ = To + + E^^ + Feoui + F,o.. 


( 12 ) 
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Table 2. Properties of nuclear matter at saturation predicted by the EoS of this work in comparison with the empirica l ranges dNewton et alj 
l2013l : [Steiner et alJ2010l : lLattimer & Liml20r3l : lHebeler et al.l2()13l : IChen et alJ2010l : lTsang et al.ll2012l : IVinas et al.ll2014l) . From left to right, the 
quantities are the energy per particle, density, incompressibility, symmetry energy, and slope parameter of the symmetry energy: L = • 

The effective nucleon mass in the present model is m* = m dBaldo et alj2008bh . The nuclear matter properties of other EoSs that will be considered 
in the sections of results are also shown. 



EjA 

(MeV) 

no 

(fm-3) 

K 

(MeV) 

^sym(^o) 

(MeV) 

L 

(MeV) 

This work 

-16.00 

0.160 

213.75 

31.92 

52.96 

Ska ILattimer & Swestv 1991) 

-16.00 

0.155 

263.18 

32.91 

74.62 

SLv4 fChabanat et al. 1998) 

-15.97 

0.160 

229.93 

32.00 

45.96 

BSk21 (Gorielv et al. 2010) 

-16.05 

0.158 

245.8 

30.0 

46.6 

TMl IShen et al. 1998a) 

-16.26 

0.145 

281.14 

36.89 

110.79 

Empirical 

-16.0 + 0.1 

0.16 + 0.01 

220 + 30 

31+2 

-60+25 


Table 3. Parameters of the Gaussian form factors and spin-orbit 
strength. 


a (fm) 

Vl (MeV) 

V[/ (MeV) 

Wo (MeV) 

0.90 

-137.024 

-117.854 

95.43 


In the calculations of open-shell nuclei we also take into account 
pairing correlations. The minimization procedure applied to the 
full functional gives a set of Hartree-like equations where the 
potential part includes in an effective way the overall exchange 
and correlation contributions. Further de t ails o f the functional 
can be found in (iBaldo et al.ll2008^ 1201 Ol l2013h . 

We note that the posited functional has only four open pa¬ 
rameters (the strengths Vl and Vu, the range a of the surface 
term, and the spin-orbit strength Wo), while the rest of the func¬ 
tional is derived from the microscopic BHF calculation. The four 
open parameters have been determined by minimizing the root- 
mean-square (rms) deviation cr^ between the theoretical and ex¬ 
perimental binding energies of a set of well-deformed nuclei in 
the rare-eart h, actinide, and supe r-heavy mass regions of the nu¬ 
clear chart (iBaldo et al.l l2008af . The optimized values of the 
parameters are given in Table [3] The predictive power of the 
functional is then tested by computing the binding energies of 
467 known spherical nuclei. A fairly reasonable rms deviation 
cte - 1.3 MeV between theory and experiment is found. In¬ 
deed, the rms deviations for binding energies and for charge 
radii of nuclei across the mass table with the BCPM functional 
are comparable to those obtained with highly successful nu¬ 
clear mean field models such as the DIS Gogny force, the SLy4 
Skyrme force, and the NL3 RMF parameter set, which for the 
same set of nuclei yield rms dev iations (Te of 1.2-1.8 MeV, see 
(IBaldo et al.ll2008al2010ll2013h . The study of ground-state de¬ 
formation properties, fission barriers, and excited octupo le states 
with the BCPM functional (iRobledo et akl I2008L1201 Ol) shows 
that the deformation properties of BCPM are similar to those 
predicted by the DIS Gogny force, which can be considered as 
a benchmark to study deformed nuclei. 

3. The outer crust 

The outer crust of a NS is the region of the star that consists of 
neutron-rich nuclei and free electrons at densities approximately 
between 10"^ g/cm^, where atoms become completely ionized, 
and 4x10^' g/cm^, where neutrons start to drip from the nuclei. 
The nuclei arrange themselves in a solid body-centreed cubic 
(bcc) lattice in order to minimize the Coulomb repulsion and are 
stabilized against /3 decay by the surrounding electron gas. At the 


low densities of the beginning of the outer crust, the Coulomb 
lattice is populated by ^®Fe nuclei. As the density of matter in¬ 
creases with increasing depth in the crust, it becomes energet¬ 
ically favourable for the system to lower the proton fraction 
through electron captures with the energy excess carried away 
by neutrinos. The system progressively evolves towards a lattice 
of more and more neutron-rich nuclei as it approaches the bot¬ 
tom of the outer crust, until the neutron drip density is reached 
and the inner crust of the NS begins. 


3. 1. Formalism for the outer crust 


To describe the structure of the outer crust we follow the for- 
malism developed by Baym, Pethick, and Suthe rland (BPS) 


(iBavmelal. 1971bh. as anplied more recently in (iRuste 

r et al.l 

120061; iRoca-Maza & Piekarewiczl 2008; Roca-Maza et al. 

2012h 

(also see (Haensel et al. 20061; Pearson et al. 2011) and refer- 


ences therein). It is considered that the matter inside the star is 
cold and charge neutral, and that it is in its absolute ground state 
in complete thermodynamic equilibrium. This is a meaningful 
assumption for non-accreting neutrons stars that have lived long 
enough to cool down and to reach equilibrium after their cre¬ 
ation. In the outer crust, the energy at average baryon number 
density n* consists of the nuclear plus electronic and lattice con¬ 
tributions: 


E{A, Z, Uh) - Ee! + Eel + El, 


(13) 


where A is the number of nucleons in the nucleus, Z is the atomic 
number, and n* - AjV (where V stands for volume). The nuclear 
contribution to Eq. (fT3l l corresponds to the mass of the nucleus, 
i.e. the rest mass energy of its neutrons and protons minus the 
nuclear binding energy: 


En(_A,Z) = M(A,Z) = (A - Z)m„ -H Zm^ - B(A,Z). (14) 


The electronic contribution reads Eei - &ei V, where the energy 
density Sei of the electrons can be considered as a that of a de¬ 
generate relativistic free Fermi gas: 


Sel 


^(2kEe + 

- ^ In [{^Fe + I Me], 


(15) 


/ 9 \l/3 

with kEe - (37r“nej being the Fermi momentum of the elec¬ 
trons, lie - {ZIA)ni, the electron number density, and nie the 
electron rest mass. The lattice energy can be written as 




(16) 
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where kfh = ^ = {AlZ^^^kfe is the average 

Fermi momentum and C; = 3.4 0665 x 10^^ for bcc lattices 
(iRoca-Maza & Piekarewiczll2008h . The lattice contribution has 
the form of the Coulomb energy in the nuclear mass formula 
with a particular prefactor and corresponds to the repulsion en¬ 
ergy among the nuclei distributed in the bcc lattice, their attrac¬ 
tion with the electrons, and the repulsion energy among the elec¬ 
trons. 

The basic assumption in the calculation is that thermal, hy¬ 
drostatic, and chemical equilibrium is reached in each layer of 
the crust. As no pressure is exerted by the nuclei, only the elec¬ 
tronic andjattice terms contribute to the pressure in the outer 
crust (iBavm et al.lll971^ . Therefore, we have 


where pe = -y Fermi energy of the electrons in¬ 

cluding their rest mass. One has to find the nucleus that at a cer¬ 
tain pressure minimizes the Gibbs free energy per particle, or 
chemical potential, p — G/A = {E — TS)IA + Pint (iBavm et al.l 
Il971hh . As far as the system can be assumed at zero tempera¬ 
ture in good approximation, since the electronic Fermi energy is 
much larger that the temperature of the star, the quantity to be 
minimized is given by 


p{A,Z,P) 


E{A,Z,nt,) P 
A lib 
M{A,Z) Z 4 
—4—+ 


( 18 ) 


For a fixed pressure, Eq. (fTSl l is minimized with respect to the 
mass number A and the atomic charge Z of the nucleus in order 
to find the optimal configuration. The nuclear masses M(A,Z) 
needed in Eq. (fTST l can be taken from experiment if they are 
available or calculated using nuclear models. 


3.2. Equation of state of the outer crust 

Although thousands of nuclear masses are experimentally de¬ 
termined not far from stability, the masses of very neutron-rich 
nuclei are not known. We used in Eq. (fTSl l the measured masses 
whenever available. Eor the unknown masses, we performed de¬ 
formed Hartree-Eock-Bogoliubov calculations with the BCPM 
energy density functional, which has been constructed from the 
microscopic BHE calculations as described in Sec. |2] We took 
the experimental data of masses from the m ost recent atomi c 
mass evaluation, i.e. the AME2012 evaluation (I Audi et al.l2012h . 
As the field of high-precision mass measurements of unstable 
neutron-rich nuclei c ontinues to advance in the radioactive beam 
facilities worldwide (S oennesse ^12013h . better constraints can 
be placed on the composition of the outer crust. Thus, we also 
included in our calculation the mass of ^^Zn recently determined 
by a Penning-trap measurement at ISOLDE-CERN (IWolf et al.l 
l2013h . Being the most neutron-rich N - 50 isotone known to 
date, ^^Zn is an impor tant nucleus in the study of the NS outer 
crust (I Wolf et al.ll2013h . 

The calculated composition of the outer crust (neutron and 
proton numbers of the equilibrium nuclei) vs., the average 
baryon density iib is displayed in Eig. [T] After the ^®Ee nu¬ 
cleus that populates the outer crust up to iib - 5 x 10“® fm“^ 
(8 X 10® g/cm^), the composition profile exhibits a sequence of 
plateaus. The effect is related with the enhanced nuclear binding 



Fig. 1. Neutron (N) and proton (Z) numbers of the predicted nuclei in 
the outer crust of a neutron star usi ng the experimental nuclear masses 
dAudi et alll2012l : rWolf et aljl2013ll when available and the BCPM en¬ 
ergy density functional or the FRDM mass formula dMbller et al.ll99^ 
for the unmeasured masses. 


for magic nucleon numbers, i.e. Z = 28 in the Ni plateau that oc¬ 
curs at intermediate densities in Eig. [T] and A = 50 and A = 82 
in the neutron plateaus that occur at higher densities. Along a 
neutron plateau, such as A = 50, the nuclei experience electron 
captures that reduce the proton number, resulting in a staircase 
structure of shells of increasingly neutron-rich isotones, until the 
jump to the next neutron plateau (A = 82) takes place. The com¬ 
position of the crust is determined by the experimental masses up 
to densities around 2.5 x 10“® fm^"* (4 x 10*° g/cm-*). At higher 
densities, model masses need to be used because the relevant nu¬ 
clei are more neutron rich and their masses are not measured. 

The results for the composition in the high-density part 
of the outer crust from the microscopic BCPM model can be 
seen in Eig. [T] They are shown along with the results from the 
macroscopic-micr oscopic finite-range droplet model (ERDM) of 
Moller, Nix et al. (iMoller et al.lll995l) that was fitted with high 
accuracy to the thousands of experimental masses available at 
the time it was published. Despite BCPM has no more than four 
fitted parameters (for reference, the sop histicated macrosco pic - 


microscopic models such as the ERDM ( 

Moller et al.lll995h and 

the Duflo-Zuker model (iDuflo & Zukeil 

1995h haye tens of ad- 


justable parameters), the predictions of BCPM are overall quite 
in consonance with those of the ERDM. In particular, the jump to 
the A = 82 plateau is predicted at practically the same density. 
Some differences are found, though, in the width of the shell 
of ^^Ni before the A = 82 plateau, and in the fact that BCPM 
reaches the neutron drip with a (very thin) shell of **"*Se, while 
the ERDM reaches the neutron drip with a shell of *'®Kr. We 
note that at the base of the outer crust (neutron drip) the baryon 
density, mass density, pressure, and electron chemical potential 
of the equilibrium configuration in BCPM are iib - 2.62 x 10“^ 
fm-^ e = 4.37 x 10“ g/cm\ P = 4.84 x IQ-"* MeV/fm^ and 
Pe = 26 .09 MeV, respectively. These valu es are close (cf. Ta¬ 
ble I of (iRoca-Maza & Piekarewiczll2008h ) to the results com- 
puted using the h ighly successful ERD M (iMoller et al.l ll99-§) 
and Duflo-Zuker (iDuflo & Zukerl[l995h models, what giyes us 
confidence that the BCPM energy density functional is well 
suited for extrapolating to the neutron-rich regions. 


Article number, page 6 of[23] 


























































B. K. Sharma et al.: Unified equation of state for neutron stars 



when available and the BCPM ene rgy density functional or the FRDM 
mass formula (iMdller et alJ 1 199.^ for the unmeasured masses. Also 
shown is the pressure from models BSk21, BPS, Lattimer-Swesty (LS- 
Ska), and Shen et al. (Shen-TMl) (see text for details). The dashed verti¬ 
cal line indicates the approximate end of the experimentally constrained 
region. 


In Fig. |2]we display the EoS, or pressure-density relation¬ 
ship, of the outer crust obtained from the experimental masses 
plus BCPM, and from the experimental masses plus the FRDM. 
One observes small jumps in the density for particular values 
of the pressure. They are associated with the change from an 
equilibrium nucleus to another in the composition. During this 
change the pressure and the chemical potential rem ain constant, 
implyi n g a finite variation of the baryon densi t y (iBavm et aT 


Ruster et'aP 120061 : iHaensel et al.l 120061 : iPearson et al 


201 11) . After the region constrained by the experimental masses 


(marked by the dashed vertical line in Fig. |2), the pressures of 
BCPM (black solid line) and the FRDM (red dashed line) ex¬ 
trapolate very similarly, with only some differences appreciated 
by the end of the outer crust. Our results for the composition and 
EoS of the outer crust are given in Table H) The quantities e and 
r in this table are, respectively, the mass density of matter and 
the adiabatic index F = (rihlP) (dPIdni,). 

Also plotted in Fig.|2]is the pressure in the outer crust from 
some popular EoSs that model the complete structure of the NS. 
The figure is drawn up to iib - 3 x 10 fm^^, thus comprising the 

change from the outer crust to the inner crust in order to allow 
comparison of the EoSs also in this region (notice, however, that 
inner-crust results are not available for the FRDM). We show in 
Fig. |2] the EoS from the recent BSk21 Skyrme nuclear effective 
force (iPearson et al.l 201^ iFantina et al.l' 201^ Potek hin et alJ 


1201iGorielv et alJ 1201 (ih tabulated in dPotekhin_e^ 


The parameters of this force were fitted (IGorielv et al 


-..EoTir 

1201 Oh to 


reproduce with high accuracy almost all known nuclear masses, 
and to various physical conditions including the neutron mat¬ 
ter EoS from microscopic calculations. We see in Fig.|2]that af¬ 
ter the experimentally constrained region, the BSk21 pressure 
is similar to the BCPM and FRDM results, with just some dis¬ 
placement around the densities where the composition changes 
from a nucleus to the next one. In the seminal work of BPS 
(iBavm et al.lll97l9) the nuclear masses for the outer crust were 


Table 4. Composition and equation of state of the outer crust. The last 
nucleus in the table with known mass is *°Zn. The experimental masses 
determine the results up to densities around 2.5 x 10“^ fm“^, while at 
higher densities the calculation of the composition involves unmeasured 
masses. 


rib 

Z 

A 

6 

P 

T 

(fm-3) 



(g cm-3) 

(erg cm“^) 


6.2203E-12 

26 

56 

1.0317E-H04 

9.5393E-H18 

1.797 

6.3129E-11 

26 

56 

1.0471E-H05 

5.3379E-H20 

1.688 

6.3046E-10 

26 

56 

1.0457E+06 

2.3241E+22 

1.586 

4.9516E-09 

26 

56 

8.2138E-(-06 

5.4155E+23 

1.470 

6.3067E-09 

28 

62 

1.0462E-H07 

7.3908E-H23 

1.459 

2.5110E-08 

28 

62 

4.1659E-H07 

5.3113E-H24 

1.400 

7.9402E-08 

28 

62 

1.3176E-(-08 

2.6112E+25 

1.369 

1.5828E-07 

28 

62 

2.6269E-t08 

6.6859E+25 

1.358 

1.6400E-07 

26 

58 

2.7220E-H08 

6.9610E-H25 

1.357 

1.7778E-07 

28 

64 

2.9508E-H08 

7.4978E-H25 

1.356 

3.1622E-07 

28 

64 

5.2496E-(-08 

1.6340E+26 

1.350 

5.0116E-07 

28 

64 

8.3212E+08 

3.0390E+26 

1.345 

7.943 lE-07 

28 

64 

1.3191E-t09 

5.6433E-)-26 

1.343 

8.5093E-07 

28 

66 

1.4132E-H09 

5.9393E-H26 

1.342 

9.2239E-07 

28 

66 

1.5319E-H09 

6.6181E-H26 

1.342 

9.9998E-07 

36 

86 

1.6609E4-09 

7.1858E-)-26 

1.341 

1.2587E-06 

36 

86 

2.0908E-(-09 

9.7829E+26 

1.340 

1.5845E-06 

36 

86 

2.6324E-H09 

1.3318E-H27 

1.340 

1.8587E-06 

36 

86 

3.0881E-H09 

1.6492E-H27 

1.339 

1.9952E-06 

34 

84 

3.3151E+09 

1.7369E+27 

1.339 

3.1620E-06 

34 

84 

5.2552E-(-09 

3.2161E+27 

1.338 

5.0116E-06 

34 

84 

8.3320E-H09 

5.9526E-H27 

1.337 

6.7858E-06 

34 

84 

1.1285E-tlO 

8.9241E-H27 

1.336 

7.5849E-06 

32 

82 

1.2615E-tlO 

9.8815E+27 

1.336 

9.9996E-06 

32 

82 

1.6635E-H10 

1.4293E+28 

1.335 

1.2589E-05 

32 

82 

2.0947E-H10 

1.9437E-H28 

1.335 

1.6595E-05 

32 

82 

2.7622E-H10 

2.8107E-H28 

1.335 

1.9053E-05 

30 

80 

3.1718E-H10 

3.2111E+28 

1.335 

2.5118E-05 

30 

80 

4.1828E-tlO 

4.6433E+28 

1.334 

3.1621E-05 

30 

80 

5.2673E-H10 

6.3132E-H28 

1.334 

3.7973E-05 

30 

80 

6.3269E-H10 

8.0596E-H28 

1.334 

4.1685E-05 

28 

78 

6.9462E+10 

8.6285E-)-28 

1.334 

5.8754E-05 

28 

78 

9.7955E-tlO 

1.3639E+29 

1.334 

6.3093E-05 

30 

84 

1.0520E-tll 

1.4867E-H29 

1.334 

7.6207E-05 

30 

84 

1.2711E-tll 

1.9125E-H29 

1.334 

8.4137E-04 

42 

124 

1.4035E-t-ll 

2.0101E+29 

1.334 

1.0964E-04 

42 

124 

1.8299E+11 

2.8616E+29 

1.334 

1.2022E-04 

40 

122 

2.0067E-t-ll 

3.1040E-H29 

1.334 

1.4125E-04 

40 

122 

2.3584E-tll 

3.8485E-H29 

1.334 

1.5667E-04 

40 

122 

2.6163E-H11 

4.4187E+29 

1.334 

1.6982E-04 

38 

120 

2.8364E-H11 

4.7062E+29 

1.334 

2.0416E-04 

38 

120 

3.4112E-H11 

6.0164E-H29 

1.334 

2.4265E-04 

38 

120 

4.0556E+11 

7.5746E-H29 

1.334 

2.6155E-04 

34 

114 

4.3721E+11 

7.7582E+29 

1.334 


provided by an early semi-empirical mass table. The correspond¬ 
ing EoS is seen to display a similar pattern with the BCPM, 
FRDM, and BSk21 results in Fig . |2] The EoS by Fattimer and 
Swestv (iF attirner & Swestvil 19911) . taken here in its Ska version 
( Fattimen2015l) (FS-Ska), and the EoS by Shen et al. (IShen et al.l 
ll998bllaHSumivoshil2015l) (Shen-TMl) were computed with, re¬ 
spectively, the Skyrme force Ska and the relativistic mean-field 
model TMl. In the two cases the calculations of masses are of 
semiclassical type and A and Z vary in a continuous way. There¬ 
fore, these EoSs do not present jumps at the densities associ¬ 
ated with a change of nucleus in the crust. Beyond this feature, 
the influence of shell effects in the EoS is rather moderate be¬ 
cause to the extent that the pressure at the densities of interest is 
largely determined by the electrons, small changes of the atomic 
number Z compared with its semiclassical estimate modify only 
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marginally the electron density and, consequently, the pressure. 
The LS-Ska EoS shows good agreement with the previously dis¬ 
cussed models, with some departure from them in the transi¬ 
tion to the inner crust. The largest discrepancies in Fig. |2] are 
observed with the Shen-TMl EoS that in this region predicts a 
softer crustal pressure with density than the other models. 


4. The inner crust: Formalism 


Below the outer crust, the inner crust starts at a density about 
4x10" g/cm^, where nuclei have become so neutron rich that 
neutrons drip in the environment, and extends until the NS core. 
The structure of the inner crust consists of a Coulomb lattice of 
nuclear clusters permeated by the gases of free neutrons and free 
electrons. This is a unique system that is not accessible to exper¬ 
iment because of the presence of the free neutron gas. The frac¬ 
tion of free neutrons increases with growing density of matter. 
At the bottom layers of the inner crust the equilibrium nuclear 
shape may change from sphere, to cylinder, slab, tube (cylindri¬ 
cal hole), and bubble (spherical hole) before going into uniform 
matter. These shapes are generically known as “nuclear pasta”. 

Full quantal Hartree-Fock (HF) calculations of the nuclear 
structures in the inner crust are complicated by the treatment of 
the neutron gas and the eventual need to deal with different ge¬ 
ometries. As a result, large scale calculations of the inner crust 
and nuclear pasta have been perform ed very often with semi- 
classical methods such as the CFDM (iFattimer & Sweswl 19911: 
iDouchin & Haenseill200lh or the Thomas-Fermi (TF) approx¬ 
imation and their varia nts, employing effective forces for the 
nucle ar interaction, see (iHaensel et al.ll20d^ IChamel & Haensel 
l2008h for reviews. Calculations of TF type, including pasta 


(Ovamatsul 19931: Gogelein & Miithed 12007 

: lOnsi et al.l 

2008 

Pearson et al.l 2012 

) and in the relativistic ( 

Cheng et al.l 

1997 

Shen et al. Il998bla 

: lAvancini et al.l 120081: iGrill et aLll2012h nu- 


clear mean-field theories. 

In this work we apply the self-consistent TF approximation 
to describe the inner crust for two reasons. First, our main aim is 
to obtain the EoS of the neutron star, which is largely driven by 
the neutron gas and where the contribution of shell corrections 
is to some extent marginal. Second, the semiclassical methods, 
as far as they do not contain shell effects, are well suited to de¬ 
scribe non-spherical shapes, i.e. pasta phases, as we shall dis¬ 
cuss below. Nevertheless, it is to be mentioned that in the case 
of spherical symmetry, shell corrections for protons in the inner 
crust may be introduced perturbatively on top of the semiclas¬ 
sical results via the Strutinsky shell-correction method. These 
corrections have been included in the cal culations of the in ¬ 
ner crust by the Bmssels-Montreal group (fC hamel et al.lj2011 
iPearson et al.l 1201 2t iFantina et al.l 1201 [ ^otekhin et al.l 120131. 
with the BSkl9-BSk21 Skyrme forces ( Gorielv et al.l 1201^ . 
Shell effects can be taken into account self-consistently using 
the HF method. Indeed, HF calculations of inner crust mat¬ 
ter were performed in spheri cal WS cells in the pion eering 
work of Negele and Vautherin (iNegele & Vautherinll 19731) . Pair¬ 
ing effects in the inner crust can have important consequences 
to describe, e.g. pulsar glitches phenomena or the cooling of 
NS. Pairing correlations have been included in BCS and HFB 
calculations of t he inner crust in the spherical WS approach , 
see for instance (iPizzochero et al. 2002t [Sandulescu et al.ll200^ 


iBaldo et^l200'^ Grill et al. 201 it iPastore et afl 201 ll) and ref¬ 

erences therein, and also in the BCS theory of anis otropic multi¬ 
band superconductivity beyond the WS approach (IChamel et al.l 
l2010h . These works are mainly devoted to investigate the super- 


fluid properties of NS inner crusts and only, to our knowledge, in 
(iBaldo et al.ll2007h the EoS of the inner crust including pairing 
correlations was reported. 

It must be mentioned as well that the formation and the 
properties of nu clear pasta have been investigated in three- 
dimensional HF ( Gogelein & Miithei 2007^ P^&_Ston3l20 1 21: 


ISchuetrumpf et al.ll2()14 ) and TF ( Okamoto et al. 2013 ) calcula¬ 

tions in cubic boxes that avoid any assumptions on the geometry 
of the system. Generally speaking, the results of these calcula¬ 
tions observe the usual pasta shapes but additional complex mor¬ 
phologies can appear a s stable or metastable state s at the transi¬ 
tions between shapes (Go gelein & Mlithei] 20071: Pais & Stone! 
1201 21: ISchuetrumpf et air 2014t lOkamoto et al.l l2013l) . Tect> 
niques based on Monte Carlo and molecular dynamics simu¬ 
lations, which do not impose a periodicity or symmetry of the 
system unlike the WS ap proximation, have also been applied to 
study nuclear pasta, se e (|HOTqwitz_eLa]J HQQ^ Watanabe et al.l 
120091: ISchneider et al.l 1201 3t H^orowitz et al.l 12015l) and reF 
erences therein. (We note that some of the quoted three- 
dimensional and simulation calculations are for pasta in super¬ 
nova matter rather than in cold neutron-star matter.) These cal¬ 
culations are highly time consuming and a detailed pressure-vs.- 
density relation for the whole inner crust is not yet available. 


4. 1. Self-consistent Thomas-Fermi description of the inner 
crust of a neutron star 

In this section we describe the method used for computing the 
structure and the EoS of the inner crust with the BCPM en¬ 
ergy density functional. Altho ugh a summary o f the approach 
has been presented elsewhere (IBaldo et al.ll2014l) . here we pro¬ 
vide a more complete report of the formalism. 

The total energy of an ensemble of A -Z neutrons, Z protons, 
and Z electrons in a spherical Wigner-Seitz (WS) cell of volume 
Vc = 4-nRll3 can be expressed as 

E — E(A,Z,Rc) - J" -I- m„n„ -I- mpHp 

-t- Eg/ {fig} -t- Seoul -1- Etgx (np,ng)jdr. (19) 

The term Ef {n„, j is the contribution of the nuclear energy 
density, without the nucleon rest masses. In our approach it reads 


E({nn,np) 


(37r2) 


2/3 


5 2m„ 

-H "V (u„(r), fipir }'), 


^^'(r) + (r) 


( 20 ) 


where the neutron and proton kinetic energy densities are writ¬ 
ten in the TF approximation and 'V {iin, fip^ is the interacting part 
provided by the BCPM nuclear energy density functional (see 
Sec.ia, which contains bulk and surface contributions. The term 
&ei {rig) in Eq. (fT9l l is the relativistic energy density due to the 
motion of the electrons, including their rest mass. Since at the 
densities prevailing in neutron-star inner crusts the Fermi energy 
of the electrons is much higher than the Coulomb energy, we 
computed £>gi using the energy density of a uniform relativistic 
Fermi gas given by Eq. (fTSl l. The term £coui {ffp, «e) in Eq. (fT^ is 
the Coulomb energy density from the direct part of the proton- 
proton, electron-electron, and proton-electron interactions. As¬ 
suming that the electrons are uniformly distributed, this term can 
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be written as 


^coul 5 

with 


Vp(r) 


f 


ehipjr') 

|r-r'| 


dr'. 


Veir) 


f 


|r-r'| 


dr'. 


a function of the cell radius Rc (or, equivalently, as a function 
of the baryon number A) and the energy differences may be be- 

(21) tween a few keV and a few eV. 

In order to obtain the EoS of the inner crust we have to com¬ 
pute the pressure. The derivation of the expression of the pres¬ 
sure in this region of the NS is given in Appe ndixiAl As shown 

(22) appendix (also see dPearson et al.l201^ '). in the inner crust 
the pressure assumes the form 


We did calculations where we allowed the profile of the elec¬ 
tron density to depend locally on position. However, we found 
this to have a marginal influence on our results for compositions 
and energies and decided to proceed with a uniform distribution 
for the electrons. Lastly, the term &exinp,ne) in Eq. (fT^ is the 
exchange part of the proton-proton and electron-electron inter¬ 
actions treated in Slater approximation: 



We consider the matter within a WS cell of radius Rc and per¬ 
form a fully variational calculation of the total energy E{A, Z, Rc) 
of Eq. ( fT9] l imposing charge neutrality and an average baryon 
density iit, - A/Vc- The fact that the nucleon densities are 
fully variational differs from some other TE calculations car- 


ried out with non-relativistic nuclear models (Ovamatsu 

1993t 

iGoeelein & Mutheill2007HOnsi et al.ll2008HPearson et al. 

I2OI2I) 

that use parametrized trial neutron and proton densities 

in the 


minimization. TE calculations of the inner crust with fully varia¬ 
tional densities using relativistic nuclear mean-field models have 


been reported in the literature (Chens etalJ 

19971 iShen et al. 

ll998hllal;lAvancini et alJl2008HC.rill et al.ll2012 

). 


Taking functional derivatives of Eq. (fT^ with respect to the 
particle densities and including the conditions of charge neutral¬ 
ity and constant average baryon density with suitable Lagrange 
multipliers, leads to the variational equations 


6'H{n„,np) 

--- -H OT„ - /i„ = 0, 

on„ 


(24) 


69{{nn,np) 

drip 


+ 

+ 


Vp(r) - U,(r) 

mp -fip^Q, 



(25) 


yjlclc + ml + y,(r) - y^(r) - j - yu, = 0, (26) 

plus the yS-equilibrium condition 


Re — Rn Rp-i (27) 

which is imposed along with the constraints of charge neutrality 
and fixed average baryon density in the cell. We note that in this 
work the chemical potentials yu„, lUp, and include the rest mass 
of the particle. 

Equations (l24li- (l27l i are solved self-consistently in a WS cell 
of spe cified size Rc following the method described in (ISil et al.l 
l2002h and the energy is calculated from Eq. (fT^ . Next, a search 
of the optimal size of the cell for the considered average baryon 
density ri/, is performed by repeating the calculation for differ¬ 
ent values of Rc, in order to find the absolute minimum of the 
energy per baryon for that n*. This can be a delicate numerical 
task because the minimum of the energy is usually rather flat as 


P = Png + Pl7" + Pll^ (28) 

where P„g is the pressure of the gas of dripped neutrons, = 

rie^ + ml - Seiirie) is the pressure of free electron gas, and 

P“ = ^ ^ corrective term from the electron 

Coulomb exchange. That is, the pressure in the inner crust is 
exerted by the neutron and electron gases in which the nuclear 
structures are embedded. On the other hand, in Appendix I bI we 
show explicitly that when the minimization of the energy per unit 
volume with respect to the radius of the WS cell is attained, sub¬ 
jected to an average density tih and charge neutrality, the Gibbs 
free energy G - PV + E per particle equals the neutron chem¬ 
ical potential, i.e. G/A = yU„. This result provides an alternative 
way to extract the pressure from the knowledge of the neutron 
chemical potential and the energy. 

The same formalism described for spherical nuclei can be ap¬ 
plied to obtain solutions for spherical bubbles (hollow spheres). 
Assuming the appropriate geometry, the method can be used in a 
similar way for other shapes (nuclear pasta), as discussed in the 
next subsection. 


4.2. Pasta phases: cylindrical and planar geometries 

The method of solving Eqs. (I24li - (l27l i can be extended to deal 
with WS cells having cylindrical symmetry (rods or “nuclear 
spaghetti”) or planar symmetry (slabs or “nuclear lasagne”). The 
length of the rods and the area of the slabs in the inner crust of 
a neutron star are effectively infinite. The calculations for these 
non-spherical geometries are simplified if one also considers the 
unit cells as rods and slabs of infinite extent in the direction per¬ 
pendicular to the base area of the rods and to the width of the 
slabs, respectively. That is, the corresponding WS cells are taken 
as rods of finite radius Rc and length L ^ oo, and as slabs of fi¬ 
nite width 2Rc and transverse area S ^ oo. With this choice of 
geometries, the number of particles and energy per unit length 
(rods) or area (slabs) are finite. By taking dV - InLrdr for rods 
and dV - Sdx for slabs in Eq. Ci, the energy calculations are 
reduced to one-variable integrals over the finite size Rc of the 
WS cell (i.e. from 0 to Rc in a circle of radius Rc in the case of 
rods, and from -Rc to +Rc along the x direction for slabs). The 
calculation of the Coulomb potentials is likewise simplified. In 
the case of rods the direct Coulomb potential can be written as 

In r I r'np(r')dr' + 

Jo 

for protons, and as 

Ve (r) = ne^ricRl ^1 - - 2 ln/?cj (30) 

for the uniform electron distribution. In the case of slabs these 
potentials read 


Vp (r) = -47re^ 




r' In r'np(r')dr' 


(29) 


Vp (x) - -4ne^ 


X I np{x')dx' + 

. Jo 




x'np(x')dx' 


(31) 
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and 

Ve (x) = -Ine^vie ^x^ + Rfj . (32) 

The other piece of the energy that depends on the geometry 
of the WS cell is the nuclear surface contribution, which is given 
by Eq. (|8]l for the BCPM energy density functional. In the case of 
spherical nuclei, performing the angular integration, the surface 
energy density of Eq. dSll can be recast as 

^surf 

X dr' 

- 7r'^^an^,(r)j. (33) 


Eor rods the nuclear surface contribution to the energy density is 
given by 


£surf [nq, Hq^ = ^ 


-tt'/V'"/' 

a 





n^ (fm‘^) 


Fig. 3. Energy per baryon of different shapes relative to uniform npe 
matter as a function of baryon density in the inner crust. 


- n^^^ariq'ir) 


(34) 


where Iq is the modified Bessel function, while for slabs it is 
given by 

^surf — —^^Q:~^Vq^qtnq{x^ 


X 


ix: 


,'{x')e 


-(x-x!}~ j a- 


dx' 


- n^^^anqi{x) 


(35) 


As happens with spherical nuclei and spherical bubbles, the 
equations corresponding to cylindrical rods can be similarly ap¬ 
plied to obtain the solutions for tube shapes (hollow rods). 


5. The inner crust: Results 

5.1. Analysis of the self-consistent solutions 

To compute the EoS of a neutron-star inner crust we need to de¬ 
termine the optimal configuration of the WS cell, i.e. the shape 
and composition that yield the minimal energy per baryon, as 
a function of the average baryon density n/,. This is done for 
each Hh value following the method explained in the previous 
section. In Eig.[2we display the results for the minimal energy 
per baryon E/A in the different shapes. The energy per baryon is 
shown relative to the value in uniform neutron-proton-electron 
(npe) matter in order to be able to appreciate the energy separa¬ 
tions between shapes. Our set of considered shapes consists of 
spherical droplets, cylindrical rods, slabs, cylindrical tubes, and 
spherical bubbles. 

As can be seen in Eig. [2 it was possible to obtain solutions 
for rods and slabs starting at densities as low as nt, ~ 0.005 fm“^. 
Solutions at lower crustal densities, i.e. from the transition den¬ 
sity between the outer crust and the inner crust until zz* = 0.005 
fm“^, were obtained for spherical droplets only. These solutions 
are not plotted in Eig. [2 for visual purposes of the rest of the 
figure. We note that at a density zz* = 0.0003 fm“^ immediately 



Hj, (fm'^) 


Fig. 4. Energy per baryon of different shapes relative to uniform npe 
matter as a function of baryon density in the high-density region of the 
inner crust. 


after the neutron drip point, i.e. at the beginning of the inner 
crust, the difference E/A - (E/A)^^ is of -2.5 MeV. Eor tube 
and bubble shapes, solutions could be obtained for crustal den¬ 
sities higher than about 0.07 fm“^, where the system is not far 
from uniform matter. 

The spherical droplet configuration is the energetically most 
favourable shape all the way up to nt ~ 0.065 fm“^, see Eig. [3] 
When the crustal density reaches ~ 0.065 fm“^ (approximately 
10'"^ g/cm^), the nucleus occupies a significant fraction of the 
cell volume and it may happen that non-sph erical structures 
have l ower energy than the spherical droplets (iRavenhall et al.l 
' l.l ll993t IOvamatsi]ll993l : iHaensel et alT 20061 : 


1 1 9831 : iLorenz et all 
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Table 5. Baryon density of the successive changes of energetically 
favourable topology in the inner crust. Units are fm“^. 


drop/rod 

rod/slab 

slab/tube 

tube/bubble 

bubble/uniform 

0.067 

0.076 

0.082 

g0.0825 

^0.0825 


IChatnel & Haenselll2008l) . We show in Fig. |4]the high-density 
region between n* = 0.07 and ni, - 0.0825 fm“^, where 
the TF-BCPM model predicts the appearance of non-spherical 
shapes as optimal configurations. Indeed, the first change of nu¬ 
clear shape occurs at iih - 0.067 fm“^ from droplets to rods. 
It can be seen in Fig.|4]that the cylindrical shape is the energeti¬ 
cally favoured configuration up to a crustal density of 0.076 fm“^ 
where a second change takes place to the planar slab shape. As 
the density of the crustal matter increases further, the energy per 
baryon of tubes and bubbles becomes progressively closer to that 
of the slabs. Very close to the crust-core transition, estimated to 
occur in the TF-BCPM model at n* = 0.0825 fm“^, there are suc¬ 
cessive shape changes from slabs to tubes at ni, = 0.082 fm“^ and 
to spherical bubbles at almost 0.0825 fm“^. At this latter density, 
the energy per baryon of the self-consistent TF bubble solution 
and that of uniform matter differ by barely -200 eV (in compar¬ 
ison, the value of EjA- m„ is of 8.3 MeV). In summary, as com¬ 
piled in Table|5j the TF-BCPM model predicts that the sequence 
of shapes in the inner crust changes from spherical droplets to 
rods, slabs, tubes, and spherical bubbles. This is in consonance 
with th e ordering of pasta phases reported in previous TF calcu - 
lations (IOvamatsulll993blAvancini et al.ll2008t iGrill et alJl2012h . 
It can be noticed that the consideration of pasta shapes renders 
the transition to the core somewhat smoother. The energy differ¬ 
ences between the most bound shape and the spherical solution 
at the same density are, though, small and do not exceed 1-1.5 
keV per nucleon. 


A typical landscape of solutions is illustrated in Fig.|5]where, 
for a density Uh — 0.077 fm“^, the energy per baryon relative to 
the neutron rest mass is shown as a function of the size of the WS 
cell for the various shapes. Usually, the curves for a given shape 
are flat around their minimum. For example, in the case of spher¬ 
ical droplets in Fig.|5j a 1% deviation of Rc from the value of the 
minimum implies a shift of merely 25 eV in EjA, but it changes 
the baryon number by as much as 25 units. Sufficiently close 
points have to be computed in order to determine precisely the Rc 
value (equivalently, the A value) that corresponds to the minimal 
energy per baryon. On the other hand, at high crustal densities 
the differences among the energy per baryon of the minima of the 
various shapes (quoted in brackets in Fig.|5]l are small; e.g. the 
minima of the slab and rod solutions in Fig. |5] differ by 190 eV. 

Figure|6]displays the cell size Rc and the proton fraction Xp - 
ZjA of the equilibrium configurations against nt, for the different 
shapes in our calculation. Rc shows a nearly monotonic down¬ 
ward trend when the density increases, in agreement with other 


studies of NS inner crusts jOyam atsul [19^ lOnsiet^ 20081 : 


iPearson et al.ll201^lAvancini et alJl2008l : lGrill et alJ 2012 ). The 
size of the cell has a significant dependence on the geometry of 
the nuclear structures, as seen by comparing Rc of the different 
shapes. In the spherical solutions, the cell radius decreases from 
almost 45 fm at rib - 0.0003 fm“^ to 13.7 fm at rib - 0.08 fm“^ 
near the transition to the core. As regards the proton fraction 
Xp, it takes quite similar values for the various cell geometries 
in the ranges of densities where we obtained solutions for the 
respective shapes. Beyond a density of the order of 0.02 fm“^ 
the proton fraction shows a weak change with density, and af¬ 



Fig. 5. The minimum energy per baryon relative to the neutron mass 
for different shapes as a function of the cell size at an average baryon 
density n* = 0.077 fm“^. The value of the absolute minimum for each 
shape is shown in MeV in brackets. 


ter rib ~ 0.05 fm ^ it smoothly tends to the value in uniform 
ripe matter. For the spherical droplet solutions, we find that Xp 
rapidly decreases from 31% at = 0.0003 fm“^ to ~3% at 
rib - 0.02 fm“^. It afterward remains almost constant, presenting 
a minimum value of 2.75% at rib — 0.045 fm“^ and then a certain 
increase up to 3.2% at the last densities before the core. 

Taking into account that except for densities close to the 
crust-core transition the spherical shape is the energetically pre¬ 
ferred configuration, and when it is not, the energy differences 
are fairly small, we restrict the discussion about the A and Z 
values in the WS cells to spherical nuclei. Figure Q depicts the 
evolution with iib of A and Z of the equilibrium spherical solu¬ 
tions. The numerical values are collected in Tabled When free 
neutrons begin to appear in the crust the number of nucleons 
contained in a cell quickly increases from A - 113 up to a max¬ 
imum A ^ 1100 at = 0.025 fm“^. From this density onwards, 
A shows a slowly decreasing trend (with a relative plateau around 
lib ~ 0.05 fm“^) till A ^ 820 at rib - 0.07 fm“^, and finally it 
presents an increase approaching the core. We see in Fig.Qthat 
the results for the atomic number Z in the inner crust show an 



n,, (fm'^) 


Fig. 6. Radius R^ of the Wigner-Seitz cell and proton fraction Xp = ZjA 
(given in percentage) in different geometries with respect to the baryon 
density. 
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n^, (fm 


Fig. 7. Mass number A (left vertical scale) and atomic number Z (right 
vertical scale) corresponding to spherical nuclei with respect to the 
baryon density. 


overall smooth decrease as a function of and a final slight 
increase before the transition to the core in agreement with the 
behaviour of A. The range of Z values lies between an upper 
value Z ^ 36 at nt - 0.01 fm“^ and a lower value Z ^ 25 
at rib - 0.07 fm“^. We note that in our calculation the mass 
and atomic numbers vary continuously and are not restricted 
to integer values because we have used the TF approximation 
which averages shell effects. The same fact explains that the 
atomic number Z of the optimal configurations does not corre¬ 
spond to proton magic numbers, as it happe ns when proton shell 
corrections are included in the calculations dNegele & Vautherinl 


1973HBaldo et al.l2007t Charnel et al.l201 ItlPe^son et al.l2012t 


Fantina et al.ll2013tlPotekhin et al.ll2013h . However, as we have 
seen from the flatness of E/A with respect to A (or Rc) around 
the minima in the discussion of Fig.|5] the resulting EoS is robust 
against the details of the composition. 

Before leaving this section, in Fig. |8] we display the spa¬ 
tial dependence of the self-consistent neutron and proton den¬ 
sity profiles for the optimal solutions in spherical WS cells with 
average baryon densities iih - 0.0475 fm“^, 0.065 fm“^, and 
0.076 fm“^. It is observed that in denser matter the size of the 
WS cell decreases, as we discussed previously, and that the 
amount of free neutrons in the gas increases, as expected. It 
can be seen that the nuclear surface is progressively washed 
out with increasing average baryon density as the nucleon dis¬ 
tributions become more uniform. At high nt the density profile 
inside the WS cell extends towards the edge of the cell, point¬ 
ing out that th e WS approximation may be close to its lim¬ 
its of validity (iNegele & Vautherinlj l973l: IChamel et^ 2007; 


Baldo et ^1200^ Pastore et al.ll201 ItlGQgelein & Miitheil 2007 


Newton & Stone l2009h . Although the proton number Z is sim¬ 

ilar for the three average baryon densities of Fig. 0 the local 
distribution of the protons is very different in the three cases. In 
Fig. [3; the proton density profile extends more than 3 fm far¬ 
ther from the origin than in Fig. [8^, while the central value of 
the proton density has decreased by more than a factor 2, hint¬ 
ing at the fact that the neutrons have a strong drag effect on the 
protons. Figure [^presents the nucleon density profiles obtained 
for cylindrical and planar geometries at the same average density 
Hb - 0.076 fm“^ as in Fig. [8};. From Figs. [8}; (droplets),!^ (rods), 
and|9}3 (slabs) we see that the size of the WS cells decreases 
with decreasing dimensionality, i.e. droplet > Rc,rod > f^c.siab- 
At high average densities near the crust-core transition, nucle¬ 


Table 6. Composition of the inner crust. 


rib 

(fm-3) 

z 

A 

Rc 

(fm) 

0.0003 

34.934 

112.991 

44.8000 

0.0005 

34.237 

153.293 

41.8300 

0.00075 

33.479 

213.369 

40.8000 

0.0010 

36.012 

264.978 

39.8450 

0.0014 

34.265 

333.809 

38.4675 

0.0017 

36.291 

376.721 

37.5400 

0.0020 

35.091 

414.026 

36.6975 

0.0025 

36.104 

466.725 

35.4550 

0.0030 

34.519 

511.212 

34.3925 

0.0035 

35.645 

550.067 

33.4775 

0.0040 

34.549 

585.320 

32.6900 

0.0050 

34.990 

648.872 

31.4075 

0.0060 

35.472 

707.137 

30.4150 

0.0075 

35.711 

787.198 

29.2625 

0.0088 

35.252 

848.825 

28.4500 

0.0100 

36.094 

898.261 

27.7825 

0.0120 

36.181 

963.496 

26.7625 

0.0135 

35.863 

999.069 

26.0450 

0.0150 

35.339 

1025.682 

25.3675 

0.0170 

34.982 

1051.388 

24.5325 

0.0180 

34.461 

1061.641 

24.1475 

0.0200 

34.036 

1078.235 

23.4350 

0.0225 

33.477 

1094.430 

22.6450 

0.0250 

32.910 

1104.446 

21.9300 

0.0275 

32.204 

1104.566 

21.2450 

0.0300 

31.290 

1092.541 

20.5625 

0.0325 

30.203 

1069.599 

19.8800 

0.0350 

29.036 

1039.295 

19.2100 

0.0375 

27.959 

1008.341 

18.5850 

0.0400 

27.152 

984.099 

18.0425 

0.0425 

26.665 

968.891 

17.5900 

0.0450 

26.549 

965.017 

17.2350 

0.0475 

26.701 

968.928 

16.9500 

0.0500 

26.955 

974.581 

16.6950 

0.0520 

27.096 

975.352 

16.4825 

0.0540 

27.065 

968.814 

16.2400 

0.0560 

26.808 

953.172 

15.9575 

0.0580 

26.322 

928.561 

15.6350 

0.0600 

25.650 

897.063 

15.2825 

0.0620 

25.080 

869.075 

14.9575 

0.0640 

24.833 

852.005 

14.7025 

0.0650 

24.750 

844.737 

14.5850 

0.0660 

24.672 

837.603 

14.4700 

0.0680 

24.613 

826.389 

14.2625 

0.0700 

24.674 

818.891 

14.0825 

0.0720 

24.875 

815.658 

13.9325 

0.0740 

25.249 

817.728 

13.8175 

0.0750 

25.502 

820.707 

13.7725 

0.0760 

25.810 

825.326 

13.7375 

0.0770 

26.190 

832.083 

13.7150 

0.0780 

26.615 

840.127 

13.7000 

0.0790 

27.111 

850.432 

13.6975 

0.0800 

27.677 

863.085 

13.7075 


ons inside the WS cell can arrange themselves in such a way 
that the region of higher density is concentrated at the edge of 
the cell, leaving the uniform region of lower density in the inner 
part of the cell. This distribution of nucleons corresponds to the 
cylindrical tube and spherical bubble configurations. In Figs.j^; 
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equal average density, the size of the cells containing tubes and 
bubbles is larger than the size of the cells accommodating rods 
and droplets, respectively, as can be appreciated by comparing 
Fig. |9^ for rods with Fig. |9}; for tubes, and Fig. [8}; for droplets 
with Fig.|9]l for bubbles. As a consequence of this fact and of the 
effectively larger value of the integration factors 2nr and Anr^ 
when the densities are accumulated near the edge of the cell, 
the total number of nucleons and the atomic number in the tube 
and bubble cells is about 1.5-2 times larger than in their rod and 
droplet counterparts. The proton fraction Xp - ZjA is, however, 
practically the same for all geometries. 


r(fm) 

Fig. 9. (a) Optimal density profile of neutrons n„ and protons n,, for 
rod shapes at average baryon density uj = 0.076 fm“^. The associated 
baryon and proton numbers, proton fraction Xp = ZjA in percentage, 
and radius of the cell are shown. The vertical dashed line marks the 
location of the end of the WS cell, (b) The same for slab shapes, (c) 
The same for tube shapes, (d) The same for bubble shapes. We note that 
the scale on the horizontal axis is the same in Figs. and|9^, and in 
Figs.|3l and[8}:. 


Il973h (label NV) and of Baldo et al. (iBaldo et al.ll2()()^ (label 


5.2. Equation of state of the inner crust 

The energy per baryon in the inner crust predicted by the BCPM 
functional is displayed against the average baryon density in 
Fig. Ho] The result is compared with other calculations available 
in the literature. This comparison will be, at the same time, use¬ 
ful to discriminate popular EoSs used in neutron-star modeling 
among each other. We show in Fig. [TOlthe results of the quan- 
tal calculations of Negele and Vautherin (iNegele & VautherinI 


Moskow). These two EoSs of the inner crust include shell effects 
and in the case of Moskow also pairing correlations. In addition 
to these models, we compare with some of the EoSs that have 
been devised to describe the comple te neutron-star stru cture. 
The EoS by Baym, Bethe, and Pethick (iBavrn et al.ll971 (la¬ 
bel B BP), the EoS by Latt imer and Swest v (iLattimer & Swestvl 
Il991h in its Ska version (lLattimei] |2015|) ( label LS-Ska), and 
the EoS by Douchin and Haensel (iDouchin & Haen^ l2001h 
(label DH-SLy4) were all obtained using the CLDM model to 
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describ e the inner crust. T he results by Shen et al. (IShen et al.l 
[199^ ISumivoshil 12015h (label Shen-TMl) were computed 
in the TF approach with trial nucleon density distributions. Fi- 
nally, in the recent EoS of the Brussels-Montreal BSk21 force 
(Pearson et alJ 201^ iFantina et al.l 1201^ iPotekhin et al.l 1201^ 
iGorielv et al.ll20ld) the inner crust was calculated in the ex¬ 
tended TF approach with trial nucleon density profiles and with 
proton shell corrections incorporated by means of the Strutinsky 
method. 

The energy in the inner crust is largely influenced by 
the properties of the neutron gas and, therefore, the EoS of 
neutron matter of the different calculations plays an ess en- 
tial role. The NV calculation (iNegele & VautherinI Il973h is 
based on a local energy density functional that closely re- 
produces the Siemens-Pandharip ande EoS of neutron matter 
(ISiemens & Pandharipande |1971[) i n the lo w-density regime. 
The Moskow calculation (IBaldo et al.l l2007h employs a semi- 
microscopic energy density functional obtained by com- 
bining the phenom enological functional of Eayans et al. 
dPayans et al.l l2000h inside the nuclear cluster with a micro¬ 
scopic part calculat ed in the Brueckner theory with the Ar- 
gonne wig potential (IWiringa et al.lll9^ to describe the neu - 
tron enyironment in the low-density r egime (|Baldo et al.ll2004l) . 
The BBP calculation (iBaym et al.l 1 1971 alio) gjyes the EoS 
based on the Br ueckner calculations for pure ne utron mat¬ 
ter o f Siemens dSiem ens & Pandharipande T97Th . The LS- 
Ska dLattimer & Swesty 19911: Lattimerl ^15h and DH-SLy4 
dPouchin & Haens^ 2001 ) EoSs were constructed using the 
Skyrme effectiye nucl ear forces Ska and SLy 4, respectiyely. The 
SLy4 Skyrme force dChabanat et al.l [l998h was parametrized, 
among other constraints, to be consistent with the microscopic 
yariational calculati on of neutron matter of Wiringa et al. 
dWiringaet al.lll988h aboye the nuclear saturation de nsity. The 
Shen-TMl EoS dShen et al.ll998blatl^umiyoshil2015h was com¬ 
puted using the relatiyistic mean field param eter set TMl for the 


19911; Ie^ 

med 

201,5|) 

Sumivoshi 

2015 

) are 


20151) are the two EoS tables in mo re widespread 


use for astrophysical simul a tions. The BSk21 EoS (Pear2on_e^^ 


2012; lEantina et al.l 1201 3^ IPotekhin et al.l 1201 3t 


Goriely et al 


20101) is based on a Skyrme force with the parameters accu¬ 


rately fitted to the known nuclear masses and constrained, among 
yarious physical conditions, to the neutron matter EoS deriyed 
within modern many-body approaches which include the contri¬ 
bution of three-body forces. 

It can be realized in Eig. [TO] that the energies per baryon 
predicted in the in ner crust by the BCPM fu nctional and by 
the NV calculation (iNegele & Vautherinlll973l) lie close oyer a 
wide range of densities, as also noti ced before (I Schuck & VinasI 
l2013l) . The result of the BBP model (iBaym et al.ll971al^ agrees 
similarly with BCPM and NV at low densities, while aboye 
rih ~ 0.02 fm“^ it predicts somewhat larger energies than BCPM 
and NV. The DH-SLy4 calculation (iDouchin & Haens'^l200ll) 
consistently predicts throughout the inner crust the largest en¬ 
ergies of all the mode ls analyzed in Pig. [TOl The energies of 


Pantina et af 


■gies ot 

JiMil 


the BSk21 calculatio n dPearscMietaD ^11- _ 

IPotekhin et al.ll20r^ iGoriely et al.ll2010t) are yery close to those 
of DH-SLy4 up to n/, ~ 0.03-0.04 fm“^. When the transition to 
the core is approached, the BSk21 energies become closer to the 
BCPM and NV results than to the DH-SLy4 result. The Moskow 
calculation (iBaldo et al.ll2007l) predicts lower energies than the 
preyious models. Howeyer, the most r emarkable differences are 
found with the results of the LS -Ska dL attimer & Swest'^ 1^91i 
lLattimerll2015l) and Shen-TMl (IShen et al.l 1998tff Sumiyoshil 



Fig. 10. Energy per baryon relative to the neutron rest mass as a func¬ 
tion of the average baryon density in the inner crust for the BCPM func¬ 
tional and other EoSs. 



Fig. 11. Pressure as a function of the average baryon density in the 
inner crust for the BCPM functional and other EoSs. The figure starts at 
= 3 X 10“'* fm“^ where Fig.|2]ended. 


12015l) calculations. It seems evident that the BCPM functional, 
as well as the results of the models constrained by some informa¬ 
tion of microscopic calculations (NV, Moskow, BBP, DI-I-SLy4, 
and BSk21), predicts overall substantially larger energies per 
baryon in the inner crust than the LS-Ska and Shen-TMl models 
that do not contain explicit information of microscopic neutron 
matter calculations. 

The pressure in the crust is an essential ingredi- 
ent entering the Tolman -Oppenheimer-Volkoff equations 
(IShaniro & Teukolskvl Il983l) that determine the mass-radius 
relation in neutron stars. The crustal pressure has also signifi¬ 
cant im plications for astrophysi cal phenomena such as pulsar 
glitches (iPiekarewicz et al.ll2014t) . As expressed in Eq. (l251 i. the 
pressure in the inner crust is provided by the free gases of the 
electrons and of the interacting dripped neutrons (aside from a 
correction from Coulomb exchange). We note, however, that 
the pressure obtained in a WS cell in the inner crust differs 
from the value in homogeneous npe matter in y6-equilibrium at 
the same average density n* owing to the lattice effects, which 
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Table 7. Equation of state of the inner crust. 


rib 

(fm-3) 

e 

(g cm-3) 

P 

(erg cm“^) 

F 

0.0003 

5.0138E-H11 

8.2141EH-29 

0.443 

0.0005 

8.3646E-H11 

1.0417EH-30 

0.560 

0.00075 

1.2555E-H12 

1.3844EH-30 

0.747 

0.0010 

1.6746E-H12 

1.6984EH-30 

0.874 

0.0014 

2.3456E-H12 

2.3837EH-30 

1.004 

0.0017 

2.8488E-H12 

2.8551EH-30 

1.070 

0.0020 

3.3522E-H12 

3.4653EH-30 

1.121 

0.0025 

4.1915E-H12 

4.4319EH-30 

1.183 

0.0030 

5.0310E-H12 

5.6159EH-30 

1.226 

0.0035 

5.8706E-H12 

6.7099EH-30 

1.257 

0.0040 

6.7106E-H12 

8.0318EH-30 

1.280 

0.0050 

8.3909E-H12 

1.0646EH-31 

1.307 

0.0060 

1.0072E-H13 

1.3476EH-31 

1.319 

0.0075 

1.2594E-H13 

1.8085EH-31 

1.322 

0.0088 

1.4781E-H13 

2.2469EH-31 

1.318 

0.0100 

1.6801E-H13 

2.6490EH-31 

1.312 

0.0120 

2.0168E-H13 

3.3595EH-31 

1.303 

0.0135 

2.2694E-H13 

3.9198EH-31 

1.299 

0.0150 

2.5221E-H13 

4.5016EH-31 

1.297 

0.0170 

2.8591E-H13 

5.2957EH-31 

1.300 

0.0180 

3.0276E-H13 

5.7163EH-31 

1.303 

0.0200 

3.3648E-H13 

6.5647EH-31 

1.314 

0.0225 

3.7864E-H13 

7.6683EH-31 

1.334 

0.0250 

4.2081E-H13 

8.8299EH-31 

1.360 

0.0275 

4.6299E-H13 

1.0060EH-32 

1.392 

0.0300 

5.0519E-H13 

1.1361EH-32 

1.427 

0.0325 

5.4740E-H13 

1.2746EH-32 

1.466 

0.0350 

5.8962E-H13 

1.4221EH-32 

1.507 

0.0375 

6.3186E-H13 

1.5794EH-32 

1.550 

0.0400 

6.7411E-H13 

1.7473EH-32 

1.594 

0.0425 

7.1637E-H13 

1.9266EH-32 

1.638 

0.0450 

7.5864E-H13 

2.1181EH-32 

1.681 

0.0475 

8.0092E-H13 

2.3227EH-32 

1.725 

0.0500 

8.4322E-H13 

2.541 IEh-32 

1.767 

0.0520 

8.7706E-H13 

2.7258EH-32 

1.800 

0.0540 

9.1092E-H13 

2.9200EH-32 

1.832 

0.0560 

9.4478E-H13 

3.1239EH-32 

1.864 

0.0580 

9.7865E-H13 

3.3370EH-32 

1.893 

0.0600 

1.0125E-H14 

3.5604EH-32 

1.922 

0.0620 

1.0464E-H14 

3.7946EH-32 

1.950 

0.0640 

1.0803E-H14 

4.0390EH-32 

1.976 

0.0650 

1.0973E-H14 

4.1651EH-32 

1.988 

0.0660 

1.1142E-H14 

4.2941EH-32 

2.000 

0.0680 

1.1481E-H14 

4.5601EH-32 

2.023 

0.0700 

1.1821E-H14 

4.8370EH-32 

2.045 

0.0720 

1.2160E-H14 

5.1245EH-32 

2.065 

0.0740 

1.2499E-H14 

5.4230EH-32 

2.083 

0.0750 

1.2669E-H14 

5.5763EH-32 

2.091 

0.0760 

1.2839E-H14 

5.7321EH-32 

2.100 

0.0770 

1.3008E-H14 

5.8903EH-32 

2.107 

0.0780 

1.3178E-H14 

6.0502EH-32 

2.114 

0.0790 

1.3348E-H14 

6.2077EH-32 

2.121 

0.0800 

1.3518E-H14 

6.3548EH-32 

2.127 


influence the electron and neutron gases. The lattice effects take 
into account the presence of nuclear structures in the crust and 
are automatically included in the self-consistent TF calculation. 

The pressure predictions in the inner crust by the BCPM 
functional are shown in Fig. [TT]in comparison with the predic¬ 


tions by the same models discussed in Fig.lTO] The initial baryon 
density in Fig. fTTIcorresponds to the last density shown in Fig. |2] 
when we studied the EoS of the outer crust in Sec. [2 In the 
inner crust, the pressure from the BCPM functional is compa¬ 
rable in general to the results of the NV, BBP, DH-SLy4, and 
BSk21 calculations. Particular agreement is observed in the in¬ 
ner crust regime between the BCPM and BSk21 pressures up 
to relatively high crustal densities. In contrast, large differences 
are found when the BCPM pressure is compared with the values 
from the Moskow, LS-Ska, and Shen-TMl models. As the crust- 
core transition is approached, these differences can be as large as 
a factor of two, which may have an influence on the predictions 
of the mass-radius relationship of neutron stars, particularly in 
small mass stars. In addition to the spherical shape, we have eval¬ 
uated the pressure for the non-spherical WS cells and the hollow 
shapes in the regime of high average baryon densities using the 
BCPM functional. However, on the one hand, as noted before, 
the pasta phases appear as the preferred configuration only in a 
relatively narrow range of densities between n* ^ 0.067 fm“^ 
and ^ 0.08 fm“^. On the other hand, the differences between 
the pressure of the spherical shape and the pressure from the 
successively favoured shapes are small, gene rally of the order o f 
1-2 keV/fm^ or less. Therefore, as we did in (iBaldo et al.ll2014l) . 
we have taken as a representative result for the whole inner crust 
the pressure calculated in the spherical droplet configurations. 
The corresponding EoS data are collected in Table Q 

An important quantity, which actually determines the re¬ 
sponse of the crust to the compression or decompression of mat¬ 
ter, is the so-called adiabatic index 


e + P dP rih dP 

P de P drib’ 


(36) 


where P is the pressure, e the mass density of matter, and Hh the 
average baryon density. In Eig. T2lwe plot the adiabatic index 
from the BCPM and DH-SLy4 dPouchin & Haenselll200lh EoS 
in the region from the last layers of the outer crust until a density 
well within the NS core (discussed in the next section). In the 
bottom layers of the outer crust the pressure is governed almost 
entirely by the ultra-relativistic electron gas, so that the value of 
F is quite close to 4/3. At the neutron drip point, F sharply de¬ 
creases by more than a factor of two due to the dripped neutrons 
that strongly soften the EoS. The just dripped neutrons contribute 
to the average baryon density but exert very little pressure. In the 
inner crust region, when the density increases the adiabatic index 
grows because the pressure of the neutron gas also increases. 
In our EoS the adiabatic index of the inner crust changes from 
F ^ 0.45 after the neutron drip point up to F ^2.1 near the 
crust-core transition (iib - 0.08 fm“^). 

In the same Fig. [T2]we report F computed in a single phase 
of homogeneous iipe matter in yS-equilibrium. It is observed that 
F in the inner crust almost coincides with the result in a sin¬ 
gle phase for densities between iib ~ 0.01 fm“^ and 0.05 fm“^. 
Above 0.05 fm“^, the adiabatic index of homogeneous npe mat¬ 
ter grows faster than in the inner crust, because the latter is soft¬ 
ened by the presence of nuclear structures and the coexistence 
between the two phases. T he predictions of the DH-SLy4 EoS 
(iDouchin & Haens'^l2001h for F in the inner crust show a simi¬ 
lar qualitative behaviour but differ quantitatively from the BCPM 
EoS. Eor example, at the bottom of the inner crust the adiabatic 
index in the DH-SLy4 EoS is F ~ 1.6, while it takes a value of 
2.1 in BCPM. This difference seems to indicate that the interact¬ 
ing part of the neutron gas at densities near the crust-core tran¬ 
sition is weaker in SLy 4 than in the BCPM functional (also see 
in this respect Fig. 1 of (iBaldo et al.ll2004h f. When the transition 
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matter is also shown. 


to the core is reached, F increases in a discontinuous way from 
2.1 to 2.5 in the BCPM EoS, due to the change from two phases 
to a single phase. With higher density in the core, the adiabatic 
index stiffens from the increase of the repulsive contributions in 
the nucleon-nucleon interaction. It is interesting to note that F 
exhibits a small sharp drop at the muon onset, which in BCPM 
is located at a density iit, - 0.13 fm“^. It arises from the ap¬ 
pearance of muons that replace some high-energy electrons and 
effectively reduce the pressure at this density. At higher densities 
F remains roughly constant, which is partly due to the increasing 
proton fraction in the npep matter of the core. 

Quasi-periodic oscillations in giant flares emitted by highly- 
magnetized neutron stars are signatures of the fundamental 


seismic shear mode, which is special ly sensitive to the nu¬ 
clear physics of neutron-star crusts dSteiner & WattsI l2009t 
ISotani et ani2012h . An important quantity for describing shear 
modes is the effective shear modulus p. It can be estimated from 
the known formula for a bcc Coulomb crystal at zero tempera- 


ture (ISteiner & Wattsll200^ISotani et al.ll20l^lStrohmaver et akl 


l1991h 


|/ = 0.1194ni^, (37) 

where Z and Rc are the proton number and the radius of a 
spherical WS cell having average baryon density nj,. In Fig. [T3] 
we display the effective shear modulus from our calculation 
with the BCPM func t ional along with the result from DH-SLy4 
dPouchin & Haens'eil 1200 ih . Because the elasticity for pasta 
phases, with the exception of spherical bubbles, is expecte d to 
be much lower than for spherical nuclei dSotani et al.ll2012l) . we 
restrict the plot in Fig.[T3]to spherical configurations, i.e. up to an 
average density ni, - 0.067 fm“^ where pasta phases start to be 
the most favourable configuration in the BCPM calculation. The 
effective shear modulus dTTl) depends on the composition of the 
crust through the proton number Z, which has a rather smooth 
variation along the inner crust with BCPM (see Fig. |7]i, and on 
the size of the WS cells that decreases from the neutron drip till 
the bottom layers of the inner crust (see Fig. The difference 
between the predictions of BCPM and SFy4 for p in Fig. [13] in¬ 
creases with density. At the higher densities of Fig. [13] the DH- 
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Fig. 13. The shear modulus from the EoS of BCPM and DH- 
SLy4 dPouchin & Haen^l200lh . 


SFy4 result is about twice the BCPM result, pointing out the 
different composition and s izes of the WS cells in the DH-SFy4 
dPouchin & Haens'dll2001h and BCPM models. Fower values of 
p for the BCPM case point in the direction of lower frequen¬ 
cies of the fundamental shear mode, but a complete analysis is 
beyond the scope of the present paper. 


6. The liquid core 

The EoS for the liquid core is deri ved in the fr amework of the 
Brueckner-Bethe-Goldstone theory dBaldoll999l) as desc ribed in 
Sec. [2 The Argonne vig potential dWiringa et al.ll 199^ is used 
as the NN interaction and three-body forces based on the so- 
called Urbana model are included in th e calculation to repro¬ 
duce the nuclear matter saturation point dSchiavilla et al.l119861: 
iBaldo et^ll997t iBaldo et^l2012tlTaranto et al.ll2013). 

In order to study the structure of the NS core, we have to cal¬ 
culate the composition and the EoS of cold, neutrino-free, cat¬ 
alyzed matter. As we discussed in the Introduction, we consider 
a NS with a core of nucleonic matter without hyperons or other 
exotic particles. We require that it contains charge neutral mat¬ 
ter consisting of neutrons, protons, and leptons (e , p~) in beta 
equilibrium, and compute the EoS for charg e neutral and beta- 
stable matter in the followi ng standard way dBaldo et al.lll997l: 
IShaniro & Teukolskvlll983]) . The Brueckner calculation yields 
the energy density of lepton/baryon matter as a function of the 
different partial densities, 

E 

e(n„. Up, He, tip) = {n„m„ + tipinp) + (n„ + np) — {n„, tip) 

+ e(np) + e(ne ), (38) 


where we have used ultrarelativistic and relativistic approx- 
imations for the energy d ensities of electrons and muons 
(IShaniro & Teukolskvll983h . respectively. In practice, it is suffi¬ 
cient to compute only the binding energy of symmetric nuclear 
matter and pure n eutron matter, since within the BHF approach it 
has been verified fealdo et alJI 199811 2000a ; ^ieune et al.ll 198^ 
IZuo et all 1 19991 : iBombaci & Fombardol 199ll) that a parabolic 
approximation for the binding energy of nuclear matter with ar¬ 
bitrary proton fraction Xp — iipliih, nt, - n„ H- n^,, is well fulfilled. 


E E 2 

— )tllj,Xp) ~ —(^Tlj^^Xp — 0.5) -t (1 '^Xp) E^ymitlh) ; 


(39) 
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Table 8. Populations of the liquid core. 


Table 9. Equation of state of the liquid core. 


lib 

(fm-3) 

Xp 

(%) 

Xe 

(%) 

Xp 

(%) 

0.0825 

2.950 

2.950 

0.000 

0.085 

3.023 

3.023 

0.000 

0.090 

3.165 

3.165 

0.000 

0.100 

3.429 

3.429 

0.000 

0.110 

3.667 

3.667 

0.000 

0.120 

3.882 

3.882 

0.000 

0.130 

4.082 

4.075 

0.007 

0.160 

4.864 

4.480 

0.384 

0.190 

5.598 

4.771 

0.827 

0.220 

6.274 

5.022 

1.253 

0.250 

6.935 

5.273 

1.662 

0.280 

7.615 

5.544 

2.071 

0.310 

8.336 

5.847 

2.489 

0.340 

9.105 

6.183 

2.922 

0.370 

9.918 

6.548 

3.371 

0.400 

10.762 

6.933 

3.829 

0.430 

11.622 

7.330 

4.292 

0.460 

12.361 

7.665 

4.697 

0.490 

13.108 

8.005 

5.102 

0.520 

13.859 

8.352 

5.507 

0.550 

14.614 

8.702 

5.912 

0.580 

15.371 

9.055 

6.315 

0.610 

16.126 

9.409 

6.717 

0.640 

16.879 

9.764 

7.115 

0.670 

17.628 

10.117 

7.510 

0.700 

18.371 

10.469 

7.902 

0.750 

19.591 

11.048 

8.542 

0.800 

20.784 

11.616 

9.167 

0.850 

21.944 

12.170 

9.775 

0.900 

23.069 

12.707 

10.363 

0.950 

24.156 

13.226 

10.930 

1.000 

25.201 

13.726 

11.475 

1.100 

27.167 

14.667 

12.501 

1.200 

28.966 

15.527 

13.439 

1.300 

30.601 

16.309 

14.292 


lib 

(fm-3) 

e 

(g cm^^) 

P 

(erg cm“^) 

E 

0.0825 

1.394EH-14 

6.916E-H32 

2.494 

0.085 

1.437EH-14 

7.454E-H32 

2.526 

0.090 

1.522EH-14 

8.626E-H32 

2.593 

0.100 

1.692EH-14 

1.138E-H33 

2.666 

0.110 

1.863EH-14 

1.472E-H33 

2.728 

0.120 

2.034EH-14 

1.869E-H33 

2.765 

0.130 

2.206EH-14 

2.333E-H33 

2.755 

0.160 

2.723E-I-14 

4.120E-H33 

2.749 

0.190 

3.244EH-14 

6.618E-H33 

2.764 

0.220 

3.770E-I-14 

9.930E-H33 

2.771 

0.250 

4.303EH-14 

1.416E-H34 

2.783 

0.280 

4.842EH-14 

1.943E-H34 

2.808 

0.310 

5.385EH-14 

2.592E-H34 

2.856 

0.340 

5.938EH-14 

3.425E-H34 

2.909 

0.370 

6.501EH-14 

4.407E-H34 

2.946 

0.400 

7.073EH-14 

5.546E-H34 

2.951 

0.430 

7.655EH-14 

6.851E-H34 

2.926 

0.460 

8.247EH-14 

8.323E-H34 

2.890 

0.490 

8.852EH-14 

9.986E-H34 

2.861 

0.520 

9.467E-I-14 

1.183E-H35 

2.840 

0.550 

1.009EH-15 

1.386E-H35 

2.821 

0.580 

1.073EH-15 

1.609E-H35 

2.803 

0.610 

1.139EH-15 

1.853E-H35 

2.786 

0.640 

1.206EH-15 

2.118E-H35 

2.771 

0.670 

1.274EH-15 

2.404E-H35 

2.756 

0.700 

1.344E-I-15 

2.712E-H35 

2.743 

0.750 

1.464EH-15 

3.275E-H35 

2.723 

0.800 

1.588EH-15 

3.902E-H35 

2.706 

0.850 

1.717EH-15 

4.595E-H35 

2.691 

0.900 

1.850EH-15 

5.357E-H35 

2.678 

0.950 

1.988EH-15 

6.190E-H35 

2.667 

1.000 

2.132EH-15 

7.095E-H35 

2.658 

1.100 

2.435EH-15 

9.135E-H35 

2.643 

1.200 

2.760EH-15 

1.149E-H36 

2.634 

1.300 

3.108EH-15 

1.418E-H36 

2.628 


where the symmetry energy iSsym can be expressed in terms of 
the difference of the energy per particle between pure neutron 
{Xp - 0) and symmetric (xp - 0.5) matter; 




1 d(E/A) 
4 dxp 


(iih, 0) 


jinb,0) - j{nh,0.5). 


(40) 


Knowing the energy density Eq. (l38l l. the various chemical 
potentials (of the species i - n, p, e,p) can be computed straight¬ 
forwardly, 



(41) 


and the equations for beta-equilibrium, 


Pi = bipn - qipe , 


(42) 


{bi and qi denoting baryon number and charge of species i) and 
charge neutrality. 


^ fiiqi = 0 , (43) 


allow one to determine the equilibrium composition n, at given 
baryon density iib and finally the EoS, 


Pirib) = nl 


d 

drib 


ejniinb)) 

lib 


de 

lib- -e = nbPn - e . 

dill, 


(44) 


In Table |8] the populations are reported for a fixed nucleon 
density, and are plotted in Eig. [14] The full red dot indicates 
the value of the nucleon density at which direct Urea processes 
set in. We remind that Ur ea processes play an important role 
in the neutron star cooling (iLattimer et al.l[l9^ lYakovlev et al.l 
l200lh . We notice that the BCPM EoS predicts a density onset 
value close to 0.53 fm“^, and therefore with our EoS medium 
mass NS can cool very quickly. In Table |9] we report the cor¬ 
responding EoS, which is represented in Eig. [TSjby a red solid 
cur ve. We notice a remarkable similarity with the EoS derived 
by (iDouchin & Haens'ell2001h (black curve), based on the ef¬ 
fective nuclear interaction SLy4 of Skyrme type. On the other 
hand, a strong difference with the Lattimer-Swesty EoS (dashed 
blue), the Shen EoS (dot-dashed, green), and the BSk21 EoS 
(dot-dashed-dashed, magenta) is observed at high densities. We 
recall that the pressures from the Lattimer-Swesty EoS and the 
Shen EoS had already been found to differ significantly from 
BCPM and SLy4 for the matter at subsaturation density in the 
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Fig. 14. The populations are displayed vs., the nucleon density for the 
BCPM EoS discussed in the text. The full red dot indicates the value of 
the nucleon density at which direct Urea processes set in. 



Fig. 15. The pressure is displayed vs., the nucleon density for 
the several EoSs discussed in the text, i.e. the BCPM (solid, red), 
the BSk21 (dot-dashed-dashed, magenta), the Lattimer-Swesty (Ska, 
dashed, blue), the Shen (dot-dashed, green), and the Douchin-Haensel 
(SLy4, dot-dot-dashed, black). A detail of the region between nt, = 
0.05 fm“^ and 0.30 fm“^ is shown in the inset. The incompressibility 
coefficients at nuclear saturation density for these models are A" = 214 
MeV (BCPM), 230 MeV (SLy4), 246 MeV (BSk21), 263 MeV (Ska), 
and 281 MeV (Shen). 


inner crust (see discussion of Fig.fTTTi. However, the BCPM and 
SLy4 pressures in the inner crust showed a concordance with 
BSk21 that remains within the core region up to about 0.2 
(see inset of Fig. [BJ but is not maintained in the extrapolation to 
higher densities, where the BSk21 and Lattimer-Swesty models 
predict the stiffest EoSs of Fig.fTSl 

Once the EoS of the nuclear m atter is known, one can solve 
the Tolman-Oppenheimer-Voikoff (IShaniro & Teukolskvlll983l) 
equations for spherically symmetric NS: 


— 

dr 

dm 

dr 




Anr~e , 


(45) 


where G is the gravitational constant, P the pressure, e the en¬ 
ergy density, m the mass enclosed within a radius r, and r the 



Fig. 16. The gravitational mass, in units of the solar mass Mq = 2 x 
lO^^g is displayed vs., the radius (left panel) and the central density 
(right panel) in units of the saturation density no = 0.16 fm“^. See text 
for details. 


(relativistic) radius coordinate. To close the equations we need 
the relation between pressure and energy density, P - P{e), 
i.e. the EoS. The integration of these equations produces the 
mass and radius of the star for given central density. It turns 
out that the mass of the NS has a maximum value as a func¬ 
tion of radius (or central density), above which the star is un¬ 
stable against collapse to a black hole. The value of the maxi¬ 
mum mass depends on the nuclear EoS, so that the observation 
of a mass higher than the maximum mass allowed by a given 
EoS simply rules out that EoS. We display in EiglT6]the rela¬ 
tion between mass and radius (left panel) and central density 
(right panel). The observed trend is consistent with the equa¬ 
tions of state displayed in EiglTS] As expected, when the stiff¬ 
ness increases the NS central density decreases for a given mass. 
The considered EoSs are compatible with the largest mass ob- 
served up to now i.e. M e, = 2.01 + 0.04 Mq in PSR J0348-I-0432 
(lAntoniadis et al.ll2013h . and displayed in Pig fTb] along with the 
previo usly observed mass of PSR J1614-2230 (iDemorestet aP 
1201 oh having Me - 1.97 + 0.04 Mq. We also notice that the 
maximum mass calculated with the BCPM and the SLy4 EoSs 
is characterized by a radius of about 10 km, which is some¬ 
what smaller than the radius calculated with the other consid¬ 
ered EoSs. Recent analyses of observations on quiescent low- 
mass X-ray bi naries (QLMXB) (Guillot & Rutled^l2014h and 
X-ray bursters (iGiiver & Ozelll^l3l) seem t o point in this direc¬ 
tion, th ough more studies could be needed (iLattimer & Steine^ 
12014bh . Eor a NS of 1.5 solar masses, the BCPM EoS predicts a 
radius of 11.63 km (see Table fTOl l. in line with the recent anal¬ 
ysis s hown in (lOzel et al.ll2015l) : see also (IChen & Piekarewi^ 
l2015h . Eor completeness, we also display in the orange hatched 
band the probability distribution for Mg and R deduced from 
five PRE (Photospheric Radius Expansion) b urst sources and five 
QLMX B sources, after a Bayesian analysis (ILattimer & Steine"^ 
l2014al) . We see that, except the Shen EoS, all EoSs are compat¬ 
ible with the observational data. High-precision determinations 
of NS radii that may be achieved in the future by planned obser¬ 
vatories su ch as the Neutron star Inte rior Composition ExploreR 
(NICER) (lArzoumanian et aT]|2014t) . should prove a powerful 
complement to maximum masses for resolving the equation of 
state of the dense matter of compact stars. 
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Table 10. Properties of the maximum mass configuration for a given 
EoS. The value of the radius R is given, as well as the radius for a star 
of mass equal to 1.5Mo. 


EoS 

^max 

R(km) 

Ri.5 (km) 

BCPM 

1.982 

9.95 

11.63 

SFy4 

2.05 

10. 

11.62 

Ska 

2.21 

10.98 

12.9 

Shen 

2.2 

12.77 

14.35 

BSk21 

2.28 

11.03 

12.57 



Fig. 17. The density profiles are shown for fixed mass configurations, 
i.e. Mg = M^ax (solid blue curve). Mg = 1.4 Mq (dashed red curve), 
and Mg = 1.0 Mq (dot-dashed green curve). Full dots indicate the onset 
of the crust. In the inset the crust thickness is displayed for each fixed 
gravitational mass. See text for details. 

In Fig. [T7] we report the density profiles for three values of 
the mass calculated with the EOS of the present work. The tran¬ 
sition density between the inner crust and the core is indicated 
by a dot along the curves. Accordingly, in the inset we report 
the ratio of the crust thickness and the total radius of the star. 
These informations can be relevant for phenomena occurring in 
the star, like glitches and deep crustal heating. 

7. Summary and outlook 

We have derived a unified equation of state for neutron stars 
with a microscopic model which is able to describe on the same 
physical framework, both the core and the crust regions. We de¬ 
scribe the neutron star structure based on modern Brueckner- 
Hartree-Fock calculations performed in symmetric and neutron 
matter. These microscopic calculations are also the basis of 
the Barcelona-Catania-Paris-Madrid energy density functional, 
devised to reproduce accurately the nuclear binding energies 
throughout the nuclear mass table. This functional is used to de¬ 
scribe the finite nuclei present in the crust of neutron stars. To 
our knowledge, this is the first time that a whole equation of 
state directly connected to microscopic results has been reported 
in the literature. 

The equation of state in the outer crust is obtained using the 
Baym-Pethick-Sutherland model, which requires the knowledge 
of atomic masses near the neutron drip line. In our calculation 
we use the experimental masses, when they are available, to¬ 
gether with the values provided by a deformed Hartree-Fock- 


Bogoliubov calculation performed with the Barcelona-Catania- 
Paris-Madrid energy density functional to estimate the unknown 
masses. We find that for average densities above e ^ 5 x 10'° 
gcm“^, where the experimental masses are unknown, the com¬ 
position of the outer crust is similar to the one predicted by the 
Finite Range Droplet Model of Moller and Nix. 

We describe the structure of the inner crust in the Wigner- 
Seitz approximation using the selfconsistent Thomas-Fermi 
method together with the Barcelona-Catania-Paris-Madrid en¬ 
ergy density functional. Electrons are considered as a relativis¬ 
tic Fermi gas with a constant density, which fill up the whole 
Wigner-Seitz cell. To obtain the optimal configuration in a cell 
for a given average density and size, the energy per baryon 
is minimized by solving self-consistently the coupled Euler- 
Fagrange equations for the neutron, proton and electron densi¬ 
ties. To obtain the most stable configuration for a given aver¬ 
age density, an additional minimization respect to the size of the 
Wigner-Seitz cell is required. 

Because the Thomas-Fermi model does not include shell cor¬ 
rections, the mass and atomic numbers corresponding to the con¬ 
figuration of minimal energy vary smoothly as a function of the 
average density. For spherical shapes, the mass number along the 
whole inner crust lie in the range between A=100 and A=800 
with a maximum around A=1100 for an average density n* ^ 
0.025 fm“^. The atomic number shows a roughly decreasing ten¬ 
dency from Z ^30 at the neutron drip up to Z ^25 at the crust 
core transition. 

Using the same Thomas-Fermi model, we have also inves¬ 
tigated the possible existence of pasta phases. To this end we 
have computed the minimal energy per baryon in Wigner-Seitz 
cells with planar and cylindrical geometries for average densities 
above ^ 0.05 fm“^ and for tube and bubble configurations 
from rib-Q-Ql fm“^ and ni=0.072 fm“^, respectively, until the 
crust-core transition density reached in our model at nh-Q.Q?i25 
fm“^. Our model predicts that up to average densities of tib = 
0.067 fm“^, spherical nuclei are the minimal energy configura¬ 
tions. With growing average densities, our model predicts the 
successive appearance of rods, slabs, tubes and spherical bub¬ 
bles as the most stable shape. 

To describe the core of neutron stars within our model we 
consider uniform matter containing neutrons, protons, electrons 
and eventually muons in yS-equilibrium, which determines the 
asymmetry of the homogeneous system. To derive the Equation 
of State in the core, we use directly the microscopic Brueckner- 
Hartree-Fock results in symmetric and neutron matter, which al¬ 
low us to easily obtain the pressure as a function of the density in 
this regime. The npe/i model is expected to be valid at least up to 
densities of about three times the saturation density, above which 
exotic matter could appear. However, we restricted ourselves to 
a NS nucleonic core and extrapolated the n pen matter to higher 
densities as was done in earlier literature dWiringa et al.lll9i^ 
iDouchin & Haens'^l2001l) . 

We have compared the predictions of our Equation of State, 
derived on a microscopic basis from the outer crust to the core, 
with the results obtained using some well known Equations of 
State available in the literature. Our Equation of State clearly 
differs from the results provided by the Fattimer and Swesty 
and Shen models, obtained in a more phenomenological way 
and widely used in astrophysical calculations. Our calculation 
agrees reasonably well in the crust with the predictions of the 
Equation of State of Douchin and Haensel and with the results 
of the BSk21 model of the Brussels-Montreal group, both based 
on Skyrme forces but including some microscopic information, 
and also shows a remarkable similarity with the former in the 
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core but differs more of the latter in this region of the neutron 
stars. 

The mass and radius of non-rotating neutron stars are ob¬ 
tained by solving the Tolman-Oppenheimer-Volkov equations. 
Our model predicts a maximal mass of about two solar masses, 
compatible with the largest mass measured up to now, and a 
radius of about 10 km. The radius obtained with our model is 
within the range of values estimated from observations of qui¬ 
escent low-mass X-ray binaries and from type I X-ray bursts. 
The mass-radius relationship computed with our model is com¬ 
parable to the results obtained using the Equation of State of 
Douchin and Haensel above the standard mass of neutron stars 
(1.4 solar masses) and differs from the predictions of the BSk21, 
Lattimer and Swesty and Shen models in this domain of neutron 
star masses. The maximal masses of the neutron stars are deter¬ 
mined by the stiffness of the Equation of State at high densities. 
As can be observed in Eigures 14 and 15, where the Equation of 
State, the maximal mass and the central density obtained using 
the different models considered in this work are displayed, when 
the stiffness increases the maximum value of the mass increases 
and, for a fixed mass, the central density decreases. 

However, there are some theoretical caveats to be consid¬ 
ered. It can be expected that quark matter appears in the centre of 
massive NS. To describe these “hybrid" NS one needs to know 
the quark matter EOS. Many models for the deconfined quark 
matter produce a too s oft EOS to support a N S of mass compati¬ 
ble with observation s jBurgio,et^ 


Maieron et al.l 120041: iNicotra et al 


l2002afl iBaldo et al.ll2b03 


I2006t iBaldo et al.l 12008; 


a 


Chen et akl 2011 ). The quark-quark interaction in the deconfined 


phase must be more repulsive in order to stiffen the EOS, and 
indeed, with a suitable quark-quark interaction, mixed quark- 
nucleon matter c an have an EOS c ompatible with two solar 
masses or more (lAlford et al.1 l2013h . An additional problem 
arises if strange matter is introduced in the NS matter. It turns 
out that BHE calculations using realistic hyperon-nucleon inter¬ 
actions known in the literature produce a too soft NS matter EOS 
and the maximu m mass is reduced to values well below the ob - 
servational limit (ISchulze et al.11200^ ISchulze & Riikenll201 ih . 
Although relativistic mean field models can be adj usted to ac¬ 
como date NS masses larger than two solar masses (lOertel et al.1 
1201 5h . and additional terms in the hyperon-nucleon i nteraction 
can provide a possible solution (lYamamoto et al.l2014l). this ”hv - 
peron puzzle" has to be considered still open (lEortin et al.l20T5h . 
In any case the nuclear EOS with the inclusion of quarks and/or 
hyperons at high density must reach a stiffness at least similar to 
the EOS with only nucleons if the two solar masses problem has 
to be overcome. 

The full BCPM EoS from the outer crust to the core in a tab¬ 
ulated form as a function of the baryon density as well as some 
other useful information is given in the text as well as supple¬ 
mentary material. There are different improvements that would 
be welcome but are left for a future work. In particular, we shall 
take into account shell effects and pairing correlations for pro¬ 
tons in order to obtain a more accurate information about the 
compositions of the different WS cells along the inner crust. This 
can be done in a perturbative way on top of the TE results by us¬ 
ing different techniques such as the Strutinsky integral method 
developed by Pearson and coworkers. On the other hand, the 
self-consistent TE-BCPM model is well suited for performing 
calculations at finite temperature, which can be of interest to in¬ 
vestigate the melting point of the crystal structure. Another nat¬ 
ural extension of this work using the TP formalism at finite tem¬ 
perature can be to obtain the EoS in the conditions of supernova 
matter. 
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Appendix A: Expression of the pressure in the 
inner crust 


In this appendix we derive the expression of the pressure in the 
model of the inner crust de scribed in Sec.ffl To this end we fol¬ 
low closely Appendix B of (iPearson et al.ll201^ . Eor the sake of 
simplicity we work in spherical symmetry, although the calcula¬ 
tion can be extended to the cylindrical and planar geometries and 
to tube and bubble configurations as well. The thermodynamical 
definition of the pressure allows one to write 


< 9 £\ _ 1 I dE 

where the volume V is identified with the volume 14 of the WS 
cell by treating the inner crust as a perfect crystal. With use of 
Leibniz’s rule for differentiation of a definite integral, the deriva¬ 
tive of the energy [see Eq. ([T9l l1 with respect to the radius Rc of 
the WS cell is 

dE 

— = 47tRI\e( (n„, + m„n„ + m,,np + Set (tie) 

"b ^coul "b ^ex ^^p-> ^ 



(A.l) 

'a,z 


+ 471 


r 


dfln dflp dfle 


r^dr. (A.2) 


Eor the calculation of the pressure this derivative is to be taken 
at constant mass and atomic numbers. If the neutron, proton, and 
electron numbers in the WS cell are fixed, we have 


0 - ^ r 

dRc Jo 

- 47rR^ fiiiRc) + 471 


4TTrni{r)dr 
■K 


Jo dRc 


(A.3) 


for i = n, p, e, which in view of Eq. (IA.2l i allows writing 
dE \ 

— = 4nRl\^ (n„, + lUntin + m,,np + Set (ric) 

"b ^coul "b &ex 


- p„n„ - ppUp - ■ (A.4) 

In the inner crust of a neutron star it is assumed that there 
are no protons in the gas of dripped neutrons and, consequently, 
lip (Rc) = 0. Therefore, at the edge of the cell we have Ef (ri„, 

= 'H{n„iRc),0) and Seoul {VpiRc) - VciRc))- 
Charge neutrality implies Seoul (Rc) = 0, because Vp (Rc) = 
Ve(Rc)- Taking into account these constraints, it is easy to 
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show that the nuclear contribution to the pressure is provided 
by only the interacting neutron gas at the edge of the WS 
cell. This neutron gas pressure is given by P„g - yu„n„(/?c) - 
[TT (nn(Rc),Q) + ninUniRc)]. The variational equation (l26l l taken 

at r - Rc implies /i^ = “ (|) ^ , and therefore 

the total pressure from the electrons is Pgi - rig yjk^g + ml - 

&ei{ne) - 5 (^) ^ We see that the free electron pressure 

P^™ - iig ^kpg + ml - Sgiitig) is modified by the electronic 

Coulomb exchange by an amount P'^J = ^ (that is, 

~ where £“ is the contribution to the energy density 

due to electronic Coulomb exchange). 

Putting together the previous results, Eq. (IA.4I) can be written 
as 

^ ^ ^ 

Thus, comparing with Eq. (lA.ll) . we see that the pressure of the 
system takes the form P - P„g + ptree + pex^ stated in Sec.|4] 


Appendix B: The minimization procedure in the 
inner crust revisited 

As alluded to in Sec.lH in this appendix we show explicitly that 
when the minimization of the energy per unit volume with re¬ 
spect to the radius Rg of the WS cell is attained, the Gibbs free 
energy per particle G/A equals the neutron chemical potential 
Pn- For simplicity we restrict the derivation to the case of the 
spherical shape. 

The optimal size Rg of the WS cell that minimizes the energy 
per unit volume under the constraints of a given average baryon 
density rih and charge neutrality is the solution of the equation 


d 

E 

1 

dV E dE 

dRg 

V 

~ V 

~ 'MgV ^ Wg 


which results in the condition 


differentiation of integrals, one has 

pRc Qy. 

A7iRl[nb{l - Xp) - n„{Rg)\ = 4;r I " r^dr, 

Jo 

7r n r*‘ dnn(r) , 

AnRllfihXp - np(Rg)] = 47r I r^dr, 

Jo 

-^r^dr. (B.4) 

0 OKg 

The expression for dEjdRg has been given in Eq. (IA.21 i. Using 
Eq. (IB.4b in Eq. (IA.2b and taking into account that np{Rg) = 0 
because there are no drip protons in the gas, one obtains 

— = 4-nRl[p{ (n„, 0 ) -l- m„n„ + Sgi {rig) “ 4 ( “ 1 

- yU„«„ -H PnHb - flgUg - llbXpijUn - Pp - ■ 

(B.5) 

Recalling that the pressure of the system is P - P„g + - 1 - P'^J, 

the result (IB.5b can be recast as 

dE 2 [ 1 

— = 47rR^ + ij„nb - ribXpi/j,, - Pp - pg)\. (B. 6 ) 

The minimization condition (IB.2b implies 

^TiRg— = 47rR^ [-P + p„nb - ribXpip,, - Pp - pg )], (B.7) 

and consequently, 

E 

P+y^Pnrib-nbXpipn-Pp-Pe)- (B. 8 ) 

When the system is in y 6 -equilibrium, the chemical potentials 
fulfill pn - Up - fig - 0, and thus Eq. (IB. 8 b gives 

PV + E^ UnA ^ G = yU„A, (B.9) 

which confirms that at the minimum of the energy per unit vol¬ 
ume with respect to the size of the cell, the Gibbs free energy 
equals the neutron chemical potential times the total number of 
baryons inside the WS cell. 
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(® 2) References 


The total number of neutrons, protons, and electrons in the WS 
cell is given by 

4;r , , 

—RgHbil - Xp) = 47r n„{r)r dr, 

^ Jo 

4;r . r^‘ 2 

—RgtibXp - 47 t ( npir)r dr, 

^ Jo 

An ^ An ^ 

—RliibXp = —Rgiig, (B.3) 

where Xp is the proton fraction. Next we take the derivative of 
(IB. 3b with respect to Rg. We note that in the present calculation, 
where we are looking for the minimum of E jV \s.. Rg, the num¬ 
ber of neutrons and protons in the WS cell is not to be assumed 
fixed, and therefore Eq. (IA.3b does not apply. Erom the deriva¬ 
tive of (IB. 3b with respect to Rg, and recalling Leibniz’s rule for 
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